Thuật toán tự khớp đường cong - P.J. Schneider

Thuật toán của Philip J. Schneider, công bố lần đầu trong cuốn Graphics Gems I (1990), là một giải pháp thích nghi (adaptive) giúp biến một tập hợp các điểm số hóa rời rạc $P_i$ thành một chuỗi các đoạn đường cong Bézier bậc ba mượt mà.

Link của bài báo tại đây.

Điểm đặc biệt của thuật toán này là việc đảm bảo tính liên tục hình học $G^1$, cho phép đường cong linh hoạt hơn so với liên tục tham số $C^1$ truyền thống.

1. Cơ sở Toán học: Đường cong Bézier bậc ba

Một đoạn Bézier bậc ba được định nghĩa bởi 4 điểm điều khiển $V_0, V_1, V_2, V_3$. Công thức tổng quát là:

Trong đó $B_i^3(t)$ là các đa thức Bernstein bậc ba:

• $B_0^3(t) = (1-t)^3$

• $B_1^3(t) = 3t(1-t)^2$

• $B_2^3(t) = 3t^2(1-t)$

• $B_3^3(t) = t^3$

Ứng với một giá trị $t$ cụ thể trong đoạn $[0,1]$ sẽ là một điểm trên đường cong.

2. Tính liên tục: $C^1$ và $G^1$ Continuity

Trong thiết kế đường cong piecewise (nối từ nhiều đoạn), việc làm thế nào để các khớp nối trông "mượt mà" là cực kỳ quan trọng. Schneider nhấn mạnh việc sử dụng $G^1$ thay vì $C^1$.

2.1. Liên tục tham số ($C^1$ - Parametric Continuity)

Một khớp nối giữa hai đoạn đường cong được gọi là $C^1$ nếu đạo hàm bậc nhất của chúng tại điểm nối là hoàn toàn bằng nhau về cả hướng và độ lớn.

• Điều kiện: $Q_1'(1) = Q_2'(0)$.

• Ý nghĩa hình học: Nếu bạn coi đường cong là quỹ đạo của một hạt đang di chuyển, $C^1$ yêu cầu hạt đó phải đi qua khớp nối với vận tốc không đổi (không tăng tốc hay giảm tốc đột ngột).

• Hạn chế: Trong khớp đường cong (curve fitting), điều này quá cứng nhắc. Nó bắt buộc các điểm điều khiển $V_1, V_2$ phải cách điểm mút một khoảng cố định bằng nhau ở cả hai phía.

2.2. Liên tục hình học ($G^1$ - Geometric Continuity)

Một khớp nối là $G^1$ nếu các vector đạo hàm bậc nhất tại điểm nối có cùng hướng nhưng có thể khác nhau về độ lớn.

• Điều kiện: $Q_1'(1) = k \cdot Q_2'(0)$ với $k > 0$.

• Ý nghĩa hình học: Các vector tiếp tuyến chỉ cần nằm trên cùng một đường thẳng (collinear). Hạt di chuyển có thể tăng hoặc giảm tốc khi đi qua khớp nối, miễn là nó không thay đổi hướng đột ngột.

• Lợi ích trong thuật toán Schneider: $G^1$ giải phóng các tham số $\alpha_1, \alpha_2$ (độ dài vector tiếp tuyến). Việc cho phép $\alpha$ thay đổi độc lập giúp đường cong có thêm "tự do" để áp sát vào các điểm dữ liệu thực tế mà vẫn đảm bảo vẻ ngoài mượt mà.

3. Các bước chính của thuật toán

Mục tiêu bài toán:

Cho tập điểm $P_0,P_1,\dots,P_n$, mục tiêu là tìm một hoặc một số các đường cong nối liền nhau sao cho các đường cong này khớp (bestfit) với tập điểm đã cho.

Bước 1: Tham số hóa ban đầu (Chord-Length Parameterization)

Để khớp đường cong, ta cần biết mỗi điểm dữ liệu $P_i$ tương ứng với giá trị $t$ nào trên đoạn $[0, 1]$. Vì chưa có đường cong, ta ước tính $t_i$ dựa trên khoảng cách giữa các điểm (độ dài dây cung hay ở đây là Chord-Length). Nghĩa là chiều dài từ điểm $P_0$ đến $P_i$, được tính bằng tổng khoảng cách các đoạn thẳng $P_0P_1,P_1P_2,\dots,P_{i-1}P_i$ và chuẩn hóa bằng cách chia cho tổng chiều dài các đoạn thẳng. Cụ thể, theo bài báo:

Sau đó chuẩn hóa $t_i = u_i / u_{n-1}$ để $t \in [0, 1]$.

Bước 2: Ràng buộc Tiếp tuyến và Alpha

Để đảm bảo tính mượt mà tại các khớp nối, hai điểm điều khiển bên trong ($V_1, V_2$) phải nằm trên hướng tiếp tuyến đơn vị $\hat{t}_1$ và $\hat{t}_2$ tại hai đầu mút ($V_0, V_3$):

• $V_1 = V_0 + \alpha_1 \hat{t}_1$

• $V_2 = V_3 + \alpha_2 \hat{t}_2$

Bài toán lúc này trở thành: Tìm $\alpha_1, \alpha_2$ để cực tiểu hóa sai số bình phương.

Bước 3: Giải Alpha bằng Phương pháp Bình phương tối thiểu (Least Squares)

Ta cần tối thiểu hóa hàm mục tiêu $S$:

Nghĩa là ở bước này, chúng ta phải tìm $\alpha_1, \alpha_2$ để hàm $S$ đạt cực tiểu. Tại điểm cực tiểu, đạo hàm riêng theo $\alpha_1, \alpha_2$ phải bằng $0$.

Đây thực chất có thể hình dung như việc di dời các điểm control points của đường cong trên trục là các vector $\hat{t_1}, \hat{t_2}$ sao cho tổng khoảng cách các điểm $P_i$ đến đường cong là nhỏ nhất.

Để giải được $\alpha$, trước tiên ta thay thế ràng buộc tiếp tuyến vào phương trình Bézier:
$V_1 = V_0 + \alpha_1 \hat{t}_1$
$V_2 = V_3 + \alpha_2 \hat{t}_2$

Thay vào $Q(t)$:

Bây giờ, ta nhóm các thành phần đã biết (hằng số) và chưa biết ($\alpha$):

Đặt:

• $f(t) = \text{ConstantPart}(t)$

• $A_1(t) = \hat{t}_1 B_1^3(t)$

• $A_2(t) = \hat{t}_2 B_2^3(t)$

Vậy: $Q(t) = f(t) + \alpha_1 A_1(t) + \alpha_2 A_2(t)$

Triển khai đạo hàm

Hàm sai số: $S = \sum_{i=0}^{n} |P_i - Q(t_i)|^2$

Đạo hàm theo $\alpha_1$:

Khai triển và chuyển vế:

Tương tự cho $\alpha_2$, ta có hệ phương trình:

Cụ thể, $X_1$ và $X_2$ đại diện cho phần sai số còn dư mà các điểm mút cố định không thể khớp được, được chiếu lên các hướng tiếp tuyến:

Ý nghĩa hình học :
$(P_i - f(t_i))$ là vector nối từ "đường cong hằng số" (đường cong khi $\alpha=0$) đến điểm dữ liệu thực tế. $X_1, X_2$ là tổng các hình chiếu của các vector sai số này lên hướng của tiếp tuyến đầu mút, trọng số bởi đa thức Bernstein.

Bước 4: Tinh chỉnh Tham số (Newton-Raphson Iteration)

Sau khi có đường cong tạm thời, các giá trị $t_i$ ban đầu (tính theo dây cung) thường không phải là vị trí tối ưu nhất. Ta cần tìm $t_i$ mới sao cho điểm $Q(t_i)$ là hình chiếu vuông góc của $P_i$ lên đường cong.

Ta giải phương trình vuông góc: $f(t) = [Q(t) - P_i] \cdot Q'(t) = 0$.
Sử dụng phương pháp Newton-Raphson:

Việc lặp lại quá trình này (khoảng 4-5 lần) giúp đường cong "hút" sát vào các điểm dữ liệu một cách đáng kể.

Bước 5: Chia tách đệ quy ($G^1$ Continuity)

Khi sai số vượt ngưỡng, ta chia đôi tại điểm lỗi nhất, tạm gọi là $P_k$. Vector tiếp tuyến chung $\vec{T}$ tại điểm chia $P_k$ được tính:

• Đoạn trái dùng $\hat{t}_2 = -\vec{T}$.

• Đoạn phải dùng $\hat{t}_1 = \vec{T}$.

Sau đó lại áp dụng fitting curve cho 2 nửa $P_0,P_1,\dots,P_k$ và $P_k,P_{k+1},\dots,P_n$

Demo