Convolution of separable filters

Đã gần hết tháng 7 và đây mới là bài đầu tiên của tháng. Bài này sẽ nói về một trường hợp đặc biệt của tích chập với các bộ lọc khả tách.

1. Convolution

Tích chập (convolution) là một thao tác quen thuộc trong xử lí ảnh. Quá trình tính toán tích chập một ảnh với một bộ lọc bất kì được minh họa trong hình sau

Minh họa các bước tính tích chập với bộ lọc Sobel.

Minh họa các bước tính tích chập với bộ lọc Sobel.

Với một bộ lọc bất kì, các thư viện xử lí ảnh hiện nay không tính tích chập một cách “ngây thơ” như trong hình, mà sử dụng biến đổi FFT để tính nhằm giảm chi phí tính toán. Tuy nhiên với một số loại bộ lọc nhất định, không nhất thiết phải dùng FFT mà vẫn có thể giảm chi phí tính toán một cách đáng kể, điển hình là các bộ lọc khả tách (separable kernels/filters) mà ta sẽ xét trong bài này.

Về mặt toán học, cho hàm f\left(x\right), tích chập của f với hàm k\left(x\right) được tính như sau:

\displaystyle r\left(x\right)=\left(f * k\right)\left(x\right)=\int f\left(x-n\right)k\left(n\right)dn

Một cách nôm na, ta “trượt” hàm k(x) lên hàm f(x) và tính tổng phần chồng lấp của 2 hàm này với nhau, do đó “hậu quả” của phép tích chập là ta được một hàm mới là sự “hoà trộn” của k(x) vào f(x). Trong trường hợp rời rạc, phép tích phân trở thành tổng:

\displaystyle r\left(i\right)=\left(f * k\right)\left(i\right)=\sum_{n} f\left(i-n\right)k\left(n\right)

Với trường hợp 2 chiều, ta chỉ cần thêm chỉ số cho chiều thứ 2:

\displaystyle r\left(i,j\right)=\left(f * k\right)\left(i,j\right)=\sum_{n}\sum_{m} f\left(i-n,j-m\right)k\left(n,m\right)

Và đây chính là cách làm việc của phép tích chập trên ảnh, được minh họa ở hình trên.

2. Separable kernels

Một bộ lọc gọi là khả tách nếu nó có thể được biểu diễn thành tích ma trận của một vector cột và một vector dòng. Chẳng hạn xét bộ lọc Sobel:

k=\left[\begin{array}{ccc} -1 & \mathbf{0} & 1\\ -2 & \mathbf{0} & 2 \\ -1 & \mathbf{0} & 1\end{array}\right] =\left[\begin{array}{c} 1 \\ 2\\ 1\end{array}\right] \times \left[\begin{array}{ccc}-1 & \mathbf{0} & 1\end{array}\right]

Như vậy bộ lọc Sobel là khả tách, do nó có thể được viết thành tích của vector cột v = \left(1, 2, 1\right)^{T} với vector dòng h = \left(-1, 0, 1\right).

Ngoài ra phép tích chập của một vector cột với một vector dòng chính là tích ma trận của 2 vector đó (dễ dàng chứng minh), nên ta có k = v \times h = v * h. Mà phép tích chập có tính kết hợp và giao hoán, nên khi tích chập một ảnh f với bộ lọc k = v*h, ta được:

\displaystyle f * k = f * \left(v * h\right) = \left(f * v\right) * h.

Nghĩa là thay vì tính tích chập theo cách cổ điển, ta có thể chia làm 2 bước: bước đầu tiên tính tích chập f * v, sau đó lấy kết quả tích chập với h.

Hiệu quả tính toán của bộ lọc khả tách

Nhưng tại sao ta lại quan tâm đến các bộ lọc khả tách? Với cách tính thông thường, một ảnh kích thước M\times N khi tích chập với một bộ lọc kích thước P\times Q, nếu không dùng FFT, thì ta sẽ cần khoảng MNPQ phép nhân và phép cộng. Nếu tách thành 2 bước như trên, ta sẽ cần khoảng MNP phép nhân và phép cộng cho bước thứ nhất, MNQ phép nhân và cộng cho bước thứ 2, tổng cộng là MN(P+Q). Như vậy bằng cách tính trong 2 bước, ta có thể tiết kiệm khoảng \frac{PQ}{\left(P+Q\right)} lần số phép tính cần thiết. Chẳng hạn với P=Q=9 thì về lí thuyết, ta có thể tăng tốc chương trình khoảng 4.5 lần.

3. Nhận biết bộ lọc khả tách

Nhưng làm thế nào để biết một bộ lọc bất kì có khả tách không? Ta biết rằng khi nhân một vector cột v với một vector dòng h, ma trận kết quả k sẽ có hạng bằng 1 (vì khi đó các cột của k là các vector phụ thuộc tuyến tính, với cơ sở là v, nghĩa là k_j = v \times h_j). Như vậy cho một bộ lọc bất kì, nếu hạng là 1 thì nó khả tách.

Theo đó, các bộ lọc trung bình, Sobel, Gaussian… đều khả tách, nhưng bộ lọc high-pass là không khả tách. Trong MATLAB có thể dễ dàng kiểm tra bằng hàm rank():


% lọc trung bình
>> a = ones(3, 3)/9;
>> rank(a)
ans =
1

% Sobel
>> b = [-1 0 1; -2 0 2; -1 0 1]
>> rank(b)
ans =
1

% Gaussian
>> c = [1 2 1; 2 4 2; 1 2 1];
>> rank(c)
ans =
1

%High-pass
>> d = [-1 -1 -1; -1 9 -1; -1 -1 -1];
>> rank(d)
ans =
2

Nhưng biết bộ lọc k là khả tách, làm thế nào để tính 2 vector vh sao cho k = v * h? Do ta đang có nhu cầu phân tích một ma trận ra thành 2 ma trận con, nên thật tự nhiên để nghĩ tới các phép phân tích ma trận: QR, SVD…

Sử dụng phân tích QR, ta có A = QR. Nếu A có hạng là p thì p cột đầu tiên của Q là độc lập tuyến tính. Do bộ lọc khả tách k có hạng bằng 1 nên nếu k = QR thì các cột của Q đều phụ thuộc tuyến tính. Hơn nữa trong ma trận tam giác R, chỉ có hàng đầu tiên có các giá trị khác 0 (Nhắc lại rằng nếu A = QR thì cột thứ j của A là tổ hợp tuyến tính của các cột trong Q, với hệ số là cột thứ j trong R: A_j=QR_j. Do đó nếu tất cả các cột trong Q đều phụ thuộc với cột đầu tiên, thì chỉ cần các hệ số trên dòng đầu tiên của R khác 0 là đủ). Do vậy ta có thể chọn v = Q\left(:, 1\right)h = R\left(1, :\right). Một ví dụ trong MATLAB sẽ cho thấy rõ điều này:


>> k = fspecial('gaussian', 5, 5)
k =
0.0369    0.0392    0.0400    0.0392    0.0369
0.0392    0.0416    0.0424    0.0416    0.0392
0.0400    0.0424    0.0433    0.0424    0.0400
0.0392    0.0416    0.0424    0.0416    0.0392
0.0369    0.0392    0.0400    0.0392    0.0369

>> [q r] = qr(k)
q =
-0.4292    0.8970   -0.0984    0.0067   -0.0386
-0.4557   -0.1715    0.4877   -0.6695   -0.2773
-0.4649   -0.3182   -0.6926    0.0772   -0.4437
-0.4557   -0.1715    0.4877    0.7238   -0.0348
-0.4292   -0.1881   -0.1869   -0.1481    0.8506
r =
-0.0859   -0.0912   -0.0931   -0.0912   -0.0859
0    0.0000    0.0000    0.0000    0.0000
0         0   -0.0000         0   -0.0000
0         0         0         0    0.0000
0         0         0         0   -0.0000

>> q(:, 1)*r(1, : )
ans =
0.0369    0.0392    0.0400    0.0392    0.0369
0.0392    0.0416    0.0424    0.0416    0.0392
0.0400    0.0424    0.0433    0.0424    0.0400
0.0392    0.0416    0.0424    0.0416    0.0392
0.0369    0.0392    0.0400    0.0392    0.0369

>> k - q(:, 1)*r(1, : )
ans =
1.0e-016 *
-0.0694   -0.2776   -0.2082   -0.2776   -0.2082
-0.0694   -0.2082   -0.2776   -0.2082   -0.2776
0   -0.2776   -0.2082   -0.2776   -0.2082
-0.0694   -0.2082   -0.2776   -0.2082   -0.2776
-0.0694   -0.3469   -0.2082   -0.3469   -0.2776

Sai số do làm tròn số thực có thể bỏ qua, ta thấy rằng k = q(:, 1)*r(1, : ).

Tất nhiên ta cũng có thể dùng phân tích SVD: k = U*S*V', với S là ma trận đường chéo chứa các trị riêng của k, và U, V là các ma trận unitary. Do k có hạng bằng 1 nên chỉ có 1 trị riêng, nghĩa là chỉ có phần tử đầu tiên của S khác o. Khi đó ta có thể chọn v = U(:, 1)*S(1, 1)h = V'(1, :)=V(:,1)'.


>> k = fspecial('gaussian', 5, 5)
k =
0.0369    0.0392    0.0400    0.0392    0.0369
0.0392    0.0416    0.0424    0.0416    0.0392
0.0400    0.0424    0.0433    0.0424    0.0400
0.0392    0.0416    0.0424    0.0416    0.0392
0.0369    0.0392    0.0400    0.0392    0.0369

>> [U S V] = svd(k)
U =
-0.4292    0.9011   -0.0614    0.0058   -0.0000
-0.4557   -0.1877    0.4020   -0.3090   -0.7071
-0.4649   -0.2757   -0.8153   -0.2076   -0.0000
-0.4557   -0.1877    0.4020   -0.3090    0.7071
-0.4292   -0.2039    0.0909    0.8752    0.0000
S =
0.2002         0         0         0         0
0    0.0000         0         0         0
0         0    0.0000         0         0
0         0         0    0.0000         0
0         0         0         0    0.0000
V =
-0.4292   -0.8690    0.2122   -0.1251         0
-0.4557    0.3297    0.4285    0.0003   -0.7071
-0.4649    0.0027   -0.4972    0.7325    0.0000
-0.4557    0.3297    0.4285    0.0003    0.7071
-0.4292    0.1658   -0.5836   -0.6691   -0.0000

>> v = U(:, 1)*S(1, 1);
v =
-0.0859
-0.0912
-0.0931
-0.0912
-0.0859

>> h = V(:, 1)'
h =
-0.4292   -0.4557   -0.4649   -0.4557   -0.4292

>> v*h
ans =
0.0369    0.0392    0.0400    0.0392    0.0369
0.0392    0.0416    0.0424    0.0416    0.0392
0.0400    0.0424    0.0433    0.0424    0.0400
0.0392    0.0416    0.0424    0.0416    0.0392
0.0369    0.0392    0.0400    0.0392    0.0369

>> k - v*h
ans =
1.0e-016 *
-0.2082   -0.3469   -0.3469   -0.4163   -0.4163
0.0694         0   -0.1388   -0.0694   -0.1388
0.0694   -0.0694   -0.0694   -0.1388   -0.1388
0.0694         0   -0.1388   -0.0694   -0.1388
0   -0.0694   -0.1388   -0.1388   -0.1388

Tất nhiên còn một cách đơn giản hơn nhiều, rất phù hợp trong các trường hợp đơn giản. Đó là nhận xét rằng k có hạng 1 nên các cột của k đều phụ thuộc tuyến tính, thành ra gọi v là cột thứ nhất của k: v = k(:, 1). Phần tử thứ i của vector h được tính bằng cách mang cột thứ i của k chia cho cột đầu tiên. Khi đó k = k(:, 1) * \left[\begin{array}{cccc} \alpha & \beta & \gamma & \ldots \end{array}\right] với alpha, beta, gamma… là kết quả các phép chia cột thứ 1, 2, 3… của k cho cột k(:, 1). Xét lại ví dụ bộ lọc Sobel ta sẽ thấy rõ điều này:

k=\left[\begin{array}{ccc} -1 & \mathbf{0} & 1\\ -2 & \mathbf{0} & 2 \\ -1 & \mathbf{0} & 1\end{array}\right] =\left[\begin{array}{c} -1 \\ -2\\ -1\end{array}\right] \times \left[\begin{array}{ccc}1 & \mathbf{0} & -1\end{array}\right]

Tất nhiên cách làm “bình dân” này chỉ phù hợp trong các trường hợp đơn giản, khi các cột của k đều có tổng tương đối bằng nhau, không có cột nào có tổng quá nhỏ để dẫn đến lỗi làm tròn số khi tính toán trên máy tính. Trong trường hợp tổng quát, một phép phân tích ma trận “stable” như QR hay SVD rõ ràng là cần thiết.

Như vậy với các bộ lọc khả tách, ta hoàn toàn có thể cài đặt một phương pháp tính tích chập tương đối hiệu quả. Hơn nữa nhận xét rằng việc tính toán tại mỗi phần từ của mảng ban đầu là hoàn toàn độc lập với nhau, nên bài toán này có thể được song song hóa một cách đơn giản. Hi vọng tôi sẽ có thời gian trình bày phiên bản song song hóa của bài toán này trên CUDA.

Advertisements

12 comments

  1. Anh V cho em hỏi là tại sao tích chập lại phải nhân xoắn như vậy ?
    Nếu được, nhờ anh giải thích giúp em tác dụng của Fourier trong việc biến miền không gian -> miền tần số ?
    Và miền tần số là gì

    Em xin cảm ơn anh 🙂

  2. Hi e,
    1. Tại sao phải nhân xoắn: theo đúng công thức thôi em 😀 Theo công thức tích chập 2 chiều, cho n và m chạy e sẽ thấy.
    2. FFT và miền tần số: Anh không nghiên cứu sâu về cái này nên chỉ nói chung chung được thôi.
    Bình thường tín hiệu được biểu diễn trong miền thời gian: cho tín hiệu X, tại thời điểm t thì X có độ lớn là h(t). Như vậy X được biểu diễn bởi hàm h(t) theo thời gian.
    Một trong những đặc tính quan trọng của tín hiệu là tần số f (tì lệ nghịch với bước sóng: f=\frac{1}{\lambda}).
    Ý tưởng của miền tần số là: Cho các tín hiệu h0(t), h1(t), h2(t).. với tần số lần lượt là f0, f1, f2… (trong Speech Recognition gọi là âm cơ bản) thì một tín hiệu h(t) bất kì đều có thể phân tích thành tổng (vô hạn) của các âm cơ bản:
    h(t) = a_0 h_0(t) + a_1 h_1(t) + a_2 h_2(t) + \cdots

    Vậy nên tín hiệu h(t) trong miền thời gian được chuyển thành dãy a_h = a_0, a_1, a_2,\cdots Dãy này là biểu diễn của h(t) trong miền tần số (với các tần số cơ bản f0, f1, f2… ở trên).

    Trong biến đổi Fourier thì các tín hiệu cơ bản h0(t), h1(t), h2(t)… đều là các hàm sin với tần số khác nhau. Như vậy biến đổi Fourier của 1 sóng là phép phân tích sóng đó thành tổng (vô hạn) của các sóng cơ bản hình sin, nhân với hệ số nào đó.
    Trong xử lí ảnh, người ta mô hình hóa tương tự, trong đó thay vì biên độ sóng h(t) theo thời gian thì ta có hàm độ sáng của ảnh theo vị trí tại mỗi pixel.

  3. Chào bạn

    Bạn cho mình hỏi: tại sao tín hiệu thực tế thì luôn luôn dương (có tín hiệu) hoặc bằng không (không có tín hiệu) nhưng khi biểu diễn dưới dạng hàm số theo thời gian f(t) và biểu diễn hình học (vẽ trên trục tọa độ) thì có cả giá trị âm?

  4. Chào bạn,
    Thực ra khi dùng kênh truyền hữu tuyến (cáp đồng, cáp quang…) thì tín hiệu trên đường truyền là tín hiệu số, chỉ là 0 hay 1. Tùy vào tín hiệu cần phát, nếu là tín hiệu analog (âm thanh, sóng truyền hình v.v..) thì sẽ cần bộ điều hợp để mã và giải mã: chuyển từ analog sang digital (ADC) ở đầu phát và từ digital sang analog (DAC) ở đầu thu.
    Tuy nhiên đây không phải chuyên môn của mình nên mình không thể nói bừa. Nếu bạn có nhu cầu thì có thể tìm thêm tài liệu về truyền tải tín hiệu số để đọc.

  5. anh Vu con o Paris ko? Dot vua roi em qua Paris choi. Tranh thu hoi bai ma may dua no hoc hoi khac em. Gio moi ve de cay mon DIP anh ah:) Anh co skype ko? Em hoi bai anh truc tiep cho de ah?

  6. Cho em hỏi là giả sử ảnh có giá trị 0 là màu đen, 255 là màu trắng, các giá trị nằm trong khoảng 0->255 là màu xám thì khi nhân tích chập, chẳng may nó có giá trị ngoài [0,255] thì sao ạ

    1. Nếu em ràng buộc giá trị trong khoảng [0, 255] thì chỉ cần clip thôi, nếu nhỏ hơn 0 thì set thành 0, nếu lớn hơn 255 thì set thành 255.

      Nhưng thông thường người ta sẽ thiết kế bộ lọc để nó không tràn.

      1. Clip là sao ạ, em đang tìm hiểu về mạng Convolution neural network và đang có thắc mắc ấy, vì tập dữ liệu rất nhiều ảnh mà filter em đoán là họ sẽ random trọng số cho nó nên chắc là thiết kế bộ lọc không tràn hơi khó ạ

      2. Ah convolutional neural net thì là chuyện khác. Các convolutional layers trong CNN không quan tâm tới giá trị của output pixel.

  7. Không quan tâm tới giá trị output của pixel là do ở softmax luôn có giá trị là nằm trong [0,1] ạ. Các giá trị trong matrận của filter là số 0 hoặc 1 ạ

    1. Thường thì giá trị trong ma trận của filter khá nhỏ, trong khoảng [-0.1, 0.1], tuỳ vào em initialize như thế nào.
      Nếu dùng convNN để classify images thì đến lớp output nó chỉ là phân phối xác suất của lớp đó vào mỗi class (nếu dùng softmax), nên không còn ý nghĩa là pixel nữa, đương nhiên không cần phải xét xem nó có phải là giá trị hợp lệ của pixel trong ảnh hay không.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s