Tính độ cong bề mặt cho phân vùng bề mặt tự do dựa trên phần mềm Matlab

Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
TRƯỜNG ĐẠI HỌC NHA TRANG • 65
THOÂNG BAÙO KHOA HOÏC
TÍNH ĐỘ CONG BỀ MẶT CHO PHÂN VÙNG BỀ MẶT TỰ DO 
DỰA TRÊN PHẦN MỀM MATLAB
SURFACE CURVATURE COMPUTATION FOR FREE-FORM SURFACE 
PARTITIONING BASED ON MATLAB PROGRAM
Nguyễn Văn Tường1 
Ngày nhận bài: 05/8/2014; Ngày phản biện thông qua: 11/8/2014; Ngày duyệt đăng: 01/12/2014
TÓM TẮT
Có thể chia bề mặt tự do thành các vùng lồi, lõm, yên ngựa khác nhau nhờ vào đặc điểm độ cong bề mặt tại các điểm 
trên bề mặt. Bài báo này trình bày việc tính độ cong của bề mặt tự do cho mục đích phân vùng. Các thông số độ cong bề 
mặt được sử dụng làm dữ liệu cho quá trình phân vùng và xác định biên các vùng khi sử dụng các đặc tính hình học của 
bề mặt và kỹ thuật mã xích trong lĩnh vực xử lý ảnh. Quá trình tính toán được thực hiện nhờ chương trình được viết bằng 
phần mềm Matlab. Dữ liệu đầu vào cho chương trình là phương trình toán học của bề mặt tự do ở dạng tường minh hoặc 
bề mặt Bspline. Tọa độ các điểm trên bề mặt cũng như trên biên của các vùng trong dữ liệu đầu ra được sử dụng cho việc 
mô hình hóa bề mặt với các vùng riêng biệt trong môi trường CAD (Computer Aided Design).
Từ khóa: Bề mặt tự do, độ cong bề mặt, phân vùng bề mặt
ABSTRACT
A free-form surface can be partitioned into different convex, concave and saddle regions thanks to the characteristics 
of surface curvatures at points on the surface. This paper presents the work of surface curvature computation for free-form 
surface partitioning. The surface curvatures are used as data for partitioning and defi ning the boundaries of regions on 
the surface when using the characteristics of surface geometry and the chain code technique in image processing fi eld. 
The computation process is performed by a Matlab program. The input data of the program are mathematical equations of 
free-form surfaces in explicit form or Bspline surfaces. The coordinates of points on the surface and on the region 
boundaries in the output data are used for modelling the surface with separate regions in CAD environment. 
Keywords: Free-form surface, surface curvatures, surface partitioning
1 TS. Nguyễn Văn Tường: Khoa Cơ khí – Trường Đại học Nha Trang
I. ĐẶT VẤN ĐỀ
Bề mặt tự do là mặt đều, trơn, thường được 
sử dụng trong các ngành thiết kế và chế tạo khuôn 
mẫu, thiết kế thân ô tô, tàu thủy, máy bay và trong 
các tác phẩm nghệ thuật. Quá trình gia công bề mặt 
tự do trên máy CNC (Computer Numerical Control) 
thường tốn nhiều thời gian do đường kính dao bị 
hạn chế bởi bán kính cong nhỏ nhất của bề mặt 
cần gia công. Một trong những phương pháp nâng 
cao năng suất gia công bề mặt tự do là chia bề mặt 
thành các vùng khác nhau và mỗi vùng có thể được 
gia công bằng các dao có đường kính khác nhau.
Chen và cs [2] đã tính các tính chất hình 
học của bề mặt tự do như độ cong Gauss, độ 
cong trung bình, độ cong cực đại và cực tiểu và
pháp véc tơ bề mặt để thành lập một véc tơ đa chiều 
cho quá trình chia vùng bề mặt tự do. Phương pháp 
nhóm cụm mờ (fuzzy clustering) và phương pháp 
C-means mờ (fuzzy C-means) đã được các tác giả 
sử dụng để chia bề mặt tự do thành các vùng lồi, 
lõm và yên ngựa. Các tính chất hình học nói trên 
cũng được Roman và cs [7, 8] để xác định biên và 
phân vùng bề mặt tự do.
Bey và cộng sự [1] đã xấp xỉ bề mặt tự do thành 
các tam giác. Các thông số pháp véc tơ và độ cong 
bề mặt được tính để xác định hình dạng cục bộ vùng 
bề mặt tại các đỉnh tam giác. Từ đó, các đỉnh này 
được nhóm thành các vùng có đồ hình khác nhau. 
Để nâng cao hiệu suất gia công bằng cách sử 
dụng dao lớn nhất có thể, Li và Zhang [4] cũng chia 
Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
66 • TRƯỜNG ĐẠI HỌC NHA TRANG
bề mặt tự do thành các vùng lồi, lõm và yên ngựa 
nhờ độ cong bề mặt. Họ gọi các vùng lõm và yên 
ngựa là những vùng tới hạn mà ở đó có thể xảy ra 
việc cắt lẹm. Các tác giả đã xây dựng chương trình 
tính toán chia vùng và kiểm tra cắt lẹm bằng ngôn 
ngữ C++.
Elber và Cohen [3] đã tiến hành phân tích độ 
cong của bề mặt tự do để nghiên cứu đặc tính đồ 
hình bề mặt. Họ đã phát triển một phương pháp lai 
sử dụng các toán tử ký hiệu và toán tử số để tính 
các độ cong bề mặt, từ đó xác định biên các đường 
biên của các vùng riêng biệt trên bề mặt. Tuy nhiên, 
các phương trình đường biên là những đa thức bậc 
cao, rất khó giải.
Nói tóm lại, cho đến nay, có nhiều công trình 
nghiên cứu phân vùng bề mặt tự do dựa trên cách 
tiếp cận dùng độ cong bề mặt. Tuy nhiên đa số các 
phương pháp đã đưa ra là khá phức tạp, khó áp 
dụng. Tường và Pokorny [9] đã đưa ra một phương 
pháp đơn giản nhưng hiệu quả để phân vùng bề mặt 
tự do. Ở đây, bề mặt tự do được chia thành các vùng 
lồi, lõm và yên ngựa dựa trên độ cong bề mặt. Đường 
biên của các vùng được xác định nhờ áp dụng kỹ 
thuật mã xích dùng trong xử lý ảnh. Quá trình tính 
toán phân vùng và xác định biên bề mặt được thực 
hiện nhờ một chương trình Matlab được viết cho bề 
mặt tự do ở dạng tường minh hoặc bề mặt Bspline. 
Bài báo này tập trung giới thiệu phương pháp tính 
toán độ cong bề mặt tự do cho mục đích phân vùng 
với sự hỗ trợ của phần mềm Matlab.
II. ĐỐI TƯỢNG NGHIÊN CỨU VÀ PHƯƠNG PHÁP 
NGHIÊN CỨU 
1. Thông số hình học của bề mặt tự do
Bề mặt tự do có thể được biểu diễn theo:
- Dạng ẩn: f(x, y, z) = 0 (1)
- Dạng tường minh: z = f(x, y) (2)
- Dạng tham số:
 S(u,v) = {Sx(u,v), Sy(u,v), Sz(u,v)} (3) 
Một số thông số hình học chủ yếu của bề mặt 
tự do là:
a. Pháp véc tơ tại một điểm
Cho một bề mặt tự do S(u, v) và một điểm bất 
kỳ trên bề mặt. Tại điểm này, Su và Sv là hai véc tơ 
tiếp tuyến theo hai phương tham số u và. Su và Sv 
không song song nhau, và một véc tơ vuông góc với 
cả hai véc tơ này là véc tơ đơn vị, được xác định bởi 
công thức (4) và được biểu diễn như trên hình 1 [6].
(4)
Véc tơ t bất kỳ vuông góc với nS được gọi là 
véc tơ tiếp tuyến với S(u,v) tại điểm p. Mặt phẳng 
chứa tất cả các véc tơ tiếp tuyến với mặt S tại điểm 
p được gọi là mặt phẳng tiếp tuyến tại điểm p, và 
được ký hiệu là Tp(S) (hình 1).
Hình 1. Pháp véc tơ đơn vị và mặt phẳng tiếp tuyến 
tại một điểm
b. Dạng toàn phương thứ nhất, F1
Dạng toàn phương thứ nhất của một bề mặt tự 
do S biểu diễn các tính chất khối của bề mặt, được 
xác định bởi [6]:
F1= dS.dS = Edu2 + 2Fdudv + Gdv2 (5)
trong đó 
 (6)
là các hệ số của dạng toàn phương thứ nhất.
c. Dạng toàn phương thứ hai, F2: 
Dạng toàn phương thứ hai mô tả độ cong của 
một bề mặt tự do, được xác định bởi [35]:
F2= −dnS .dS = Ldu2 + 2Mdudv + Ndv2 (7)
trong đó 
 (8)
là các hệ số của dạng toàn phương thứ hai.
d. Độ cong Gauss (K) và độ cong trung bình (H)
Cho một bề mặt tự do S(u,v) và p là một điểm 
bất kỳ trên nó. Gọi (Q) là mặt phẳng chứa pháp véc 
tơ của mặt S tại điểm p. Giao tuyến của Q là S là 
một đường cong có độ cong nhất định (hình 2). Khi 
mặt (Q) quay xung quanh pháp véc tơ nói trên thì 
độ cong của đường cong thay đổi. Ơ-le đã chính 
minh rằng tồn tại các hướng mà ở đó độ cong của 
đường cong đạt đạt giá trị cực tiểu và cực đại [6]. 
Các độ cong ở các hướng này được gọi là các độ 
cong chính tắc và các hướng độ cong chính tắc 
vuông góc nhau. 
Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
TRƯỜNG ĐẠI HỌC NHA TRANG • 67
Các điểm không có tính chất này được mã hóa 
thành số 0.
- Nếu K ≥ 0 và H > 0: lưu điểm vào ma trận các 
điểm vùng lõm, mã hóa điểm lưu thành số 2. Các 
điểm không có tính chất này được mã hóa thành số 0.
- Nếu K < 0: lưu điểm vào ma trận các điểm 
vùng yên ngựa, mã hóa điểm lưu thành số 3. Các 
điểm không có tính chất này được mã hóa thành 
số 0.
Cấu trúc ma trận các vùng riêng biệt nêu trên 
tương tự như cấu trúc của ma trận của ảnh nhị 
phân. Điều này tạo điều kiện dễ dàng cho việc áp 
dụng kỹ thuật mã xích trong xử lý ảnh để xác định 
biên các vùng đã phân. Các điểm biên sẽ được 
dùng để xây dựng các đường cong không gian ba 
chiều dùng cho việc chia bề mặt tự do thành các 
vùng khác nhau trong môi trường CAD.
Để thực hiện tính toán theo thuật toán nói trên 
và căn cứ các phương trình (9) đến (12), các bước 
để tính độ cong của một bề mặt tự do được thực 
hiện như sau:
- Tạo tập hợp các điểm trên bề mặt tự do 
đã cho.
- Tính toán các pháp véc tơ đơn vị tại tất cả các 
điểm trên bề mặt. 
- Tính các hệ số dạng toàn phương thứ nhất 
và thứ hai.
- Tính độ cong Gaussian, độ cong trung bình và 
các độ cong chính tắc.
Trong nghiên cứu này, việc tính toán phân vùng 
và xác định biên các vùng được thực hiện bằng một 
chương trình Matlab. Chương trình gồm các tập tin 
M-function và M-script để tạo mô hình toán của bề 
mặt, tính độ cong bề mặt, phân vùng và xác định 
biên các vùng. Hàm tính toán độ cong có nội dung 
cơ bản như sau:
- Tính đạo hàm bậc nhất và đạo hàm bậc hai 
theo các biến u và v bằng cách sử dụng hàm tiêu 
chuẩn gradient.
[Xu,Xv] = gradient(X); 
[Yu,Yv] = gradient(Y);
[Zu,Zv] = gradient(Z); 
[Xuu,Xuv] = gradient(Xu); 
[Yuu,Yuv] = gradient(Yu);
[Zuu,Zuv]= gradient(Zu);
[Xuv,Xvv] = gradient(Xv);
[Yuv,Yvv] = gradient(Yv);
[Zuv,Zvv] = gradient(Zv);
X, Y và Z ở đây là những mảng 2 chiều của 
các điểm trên bề mặt. Những mảng này phải được 
chuyển thành các véc tơ để thực hiện các phép tính 
về sau.
Hình 2. Độ cong của bề mặt tự do
Độ cong Gaussian (K), độ cong trung bình (H) và 
các độ cong chính tắc Kmin và Kmax của bề mặt S(u, v), 
tại điểm p, được tính bằng các công thức [6]:
(9)
 (10)
 (11)
 (12)
Dựa vào các giá trị của độ cong Gauss, độ cong 
trung bình và các độ cong chính tắc, các điểm trên 
một bề mặt tự do có thể được chia thành sáu loại 
như sau [5]:
- Điểm eliptic lõm: Nếu K > 0 và H > 0.
- Điểm eliptic lồi: Nếu K > 0 và H < 0.
- Điểm hyperbolic: Nếu K <0.
- Điểm parabolic lõm: Nếu K = 0 và H > 0.
- Điểm parabolic lồi: Nếu K = 0 và H < 0.
- Điểm rốn phẳng: Nếu K = 0 và H = 0.
Để chia một bề mặt tự do thành các vùng lồi (kể 
cả vùng phẳng), vùng lõm và vùng yên ngựa, hình 
dạng bề mặt cục bộ quanh một điểm có thể được 
chia thành ba loại vùng khác nhau như sau [5]:
* K ≥ 0 và H £ 0: hình dạng bề mặt cục bộ lồi.
* K ≥ 0 và H > 0: hình dạng bề mặt cục bộ lõm.
* K < 0 và H ¹ 0: hình dạng bề mặt cục bộ yên ngựa.
2. Tính độ cong của bề mặt tự do
Trong nghiên cứu này, thuật toán phân vùng bề 
mặt tự do như sau:
(a) Tạo tập hợp lưới điểm bề mặt {p} từ mô hình 
toán học của bề mặt S và lưu tất cả các điểm vào
một ma trận chung.
(b) Tính các thông số K and H tại mỗi điểm pi,j.
(c) Xét mỗi điểm pi,j thuộc tập {p}:
- Nếu K ≥ 0 và H £ 0: lưu điểm vào ma trận 
các điểm vùng lồi, mã hóa điểm lưu thành số 1. 
Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
68 • TRƯỜNG ĐẠI HỌC NHA TRANG
III. KẾT QUẢ NGHIÊN CỨU VÀ THẢO LUẬN
Trong nghiên cứu này, chương trình Matlab 
được viết cho tường minh hoặc bề mặt Bspline. Bài 
báo này sử dụng bề mặt tự do ở dạng tường minh 
để minh họa cho việc tính toán độ cong bề mặt.
Ví dụ: Cho bề mặt tự do được biểu diễn bởi 
phương trình , trong đó x và y có 
giá trị trong đoạn [-1,3]. Giả sử ma trận điểm lưới 
cần tạo trên bề mặt có cỡ là 41´41 theo hai phương 
x và y. 
Trong nghiên cứu này, chương trình Matlab 
được chạy trên máy tính xách tay (Intel Core i5, 
1,80GHz, RAM 4 GB) cài đặt hệ điều hành Windows 
7. Hình 3 và hình 4 trình bày kết quả tính K, H, Kmax 
và Kmin tại một số điểm lưới trên bề mặt.
- Tính pháp véc tơ tại các điểm trên bề mặt:
 m = cross(Xu,Xv,2); q = sqrt(dot(m,m,2)); n = m./[q q q];
- Tính các hệ số dạng toàn phương thứ nhất:
 E = dot(Xu,Xu,2); F = dot(Xu,Xv,2); G = dot(Xv,Xv,2);
- Tính các hệ số dạng toàn phương thứ hai:
 L = dot(Xuu,n,2); M = dot(Xuv,n,2); N = dot(Xvv,n,2);
- Tính độ cong Gauss:
 K = (L.*N - M.^2)./(E.*G - F.^2);
- Tính độ cong trung bình:
 H = (E.*N + G.*L - 2.*F.*M)./(2*(E.*G - F.^2));
- Tính độ cong chính tắc:
 Kmax = H + sqrt(H.^2 - K); Kmin = H - sqrt(H.^2 - K);
Hình 3. Minh họa giá trị K, H tại một số điểm lưới
Hình 4. Minh họa giá trị Kmax, Kmin tại một số điểm lưới
Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
TRƯỜNG ĐẠI HỌC NHA TRANG • 69
Kết quả tính toán cho thấy bề mặt đã cho có 
2 vùng lõm, 1 vùng lồi và 3 vùng yên ngựa. Chương 
trình cũng cho kết quả tọa độ các điểm biên của 
6 vùng này. Tuy nhiên, đối với ví dụ này, chỉ cần sử 
dụng các điểm biên của các vùng lõm và vùng lồi để 
xây dựng hai đường cong không gian cho mục đích 
chia vùng bề mặt trong môi trường CAD. Trên hình 5
là mô hình bề mặt đã cho với các điểm biên của 
các vùng lõm và vùng lồi, được hiển thị trong môi 
trường Matlab.
Hình 5. Mô hình bề mặt với các điểm biên của vùng lồi, lõm
Các điểm biên của vùng lồi được chứa trong 
ma trận PB24, gồm 71 điểm (hình 6a). Ma trận 
PB13cell{1,1} trên hình 6b là ma trận chứa 27 điểm 
biên của vùng lõm thứ nhất và ma trận PB13cell{2,1} 
trên hình 6c chứa 60 điểm biên của vùng lõm thứ 
hai. Trong các ma trận này, các cột 1, 2 và 3 tương 
ứng với tọa độ x, y và z của các điểm. 
Trong chương trình tính toán này, các giá trị 
tọa độ của các điểm trên bề mặt (ở một ma trận 
riêng) cũng như trên các biên có thể dễ dàng được 
lưu ở dạng fi le Excel. Điều này tạo thuận lợi cho 
việc nhập dữ liệu xây dựng bề mặt đã cho và các 
đường cong biểu diễn biên của các vùng trên bề 
mặt bằng các phần mềm CAD như Catia, Creo, NX,
SolidWorks, Trong môi trường CAD, các đường 
cong biên sẽ được sử dụng làm công cụ xén bề mặt 
để chia bề mặt nguyên thành các vùng riêng biệt. 
Để tăng độ chính xác của các đường cong biên, 
có thể tạo ma trận điểm lưới bề mặt với mức độ 
điểm dày hơn. Khi đó số lượng các điểm trên các 
biên sẽ tăng. Tuy nhiên, lúc này máy tính sẽ tính 
toán lâu hơn. Ở ví dụ này, với ma trận điểm 41´41 
thì thời gian chạy chương trình là 0,35 giây. Thời 
gian chạy chương trình khi ma trận điểm 201´201 
là 5,7 giây. Bề mặt có kích thước càng lớn thì thời 
gian chạy chương trình của máy tính sẽ càng lâu.
 (a) (b) (c)
Hình 6. Minh họa tọa độ các điểm biên trên các vùng lồi (a) và lõm (b, c)
Tạp chí Khoa học - Công nghệ Thủy sản Số 4/2014
70 • TRƯỜNG ĐẠI HỌC NHA TRANG
Chương trình Matlab trong nghiên cứu này mới 
chỉ thực hiện tính toán phân tất cả các vùng hiện 
có trên bề mặt đã cho thành các vùng riêng biệt. 
Chương trình chưa thực hiện việc tối ưu hóa quá 
trình phân vùng như kết hợp các vùng có diện tích 
quá nhỏ, có bán kính cong nhỏ nhất khác biệt không 
đáng kể so với bán kính cong nhỏ nhất của các 
vùng liền kề, thành một vùng lớn hơn nhằm nâng 
cao năng suất gia công. Hiện tại, việc kết hợp các 
vùng này có thể được thực hiện thủ công trong môi 
trường CAD bằng cách không sử dụng biên các 
vùng có diện tích nhỏ để trong quá trình chia mặt. 
IV. KẾT LUẬN VÀ KIẾN NGHỊ
Trong nghiên cứu này, các thông số độ cong 
của bề mặt được tính toán để phục vụ cho việc 
phân vùng bề mặt tự do thành các vùng lồi, lõm 
và yên ngựa. Việc tính toán được thực hiện bằng 
chương trình Matlab. Dữ liệu đầu vào trong nghiên 
cứu này là phương trình toán học của bề mặt tự 
do dạng tường minh hoặc bề mặt Bspline. Dữ liệu 
đầu ra là tập hợp các điểm trên bề mặt và các điểm 
trên biên các vùng của bề mặt này, để phục vụ cho 
việc xây dựng bề mặt tự do với các vùng riêng biệt 
trong môi trường CAD. Kết quả cho thấy chương 
trình đã tạo được các dữ liệu cần thiết cho việc phân 
vùng bề mặt phức tạp. Tuy nhiên, để nghiên cứu 
này hoàn thiện hơn, cần thực hiện tối ưu hóa quá 
trình phân vùng theo tiêu chí thời gian gia công bé 
nhất trong trường hợp đã chọn đúng các dao khác 
nhau để gia công các vùng.
TÀI LIỆU THAM KHẢO
1. Bey M., Bendifallah M., Kader S.and Boukhalfa K., 2008. Cutting tool combination and machining strategy affectation based 
on the determination of local shapes for free form surfaces. International Conference on Smart Manufacturing Application, 
120-125.
2. Chen Z.C., Dong Z. and Vickers G.W., 2003. Automated surface subdivision and tool path generation for 3½½-axis CNC 
machining of sculptured parts. Computers in Industry, 50(3), 319-331.
3. Elber G. and Cohen E., 1993. Second order surface analysis using hybrid symbolic and numeric operators. Transactions on 
Graphics, 12(12), 60-178.
4. Li L. L. and Zhang Y. F., 2004. Cutter selection for 5-axis milling based on surface decomposition. 8th International 
Conference on Control, Automation, Robotics and Vision, Kunming, China, 3, 1863- 1868. 
5. Makhanov S.S. and Anotaipaiboon W., 2007. Advanced numerical methods to optimize cutting operations of fi ve-axis milling 
machines. Springer-Verlag, Berlin.
6. Radzevich S..P., 2008. CAD/CAM of Sculptured surfaces on multi-axis NC machine: The DG/K-based approach. Morgan 
& Claypool, USA.
7. Roman A, Bedi S. and Ismail F. Three-half and half-axis patch-by-patch NC machining of sculptured surfaces. International 
Journal of Advanced Manufacturing Technology, 2006, Vol 29, No 5-6, 524-531.
8. Roman A. Surface partitioning for 3+2-axis Machining. D.Phil. Thesis, University of Waterloo, Canada, 2007.
9. Tuong N.V. and Pokorny P., 2010. A practical approach for parttioning free-form surfaces. International Journal of Computer 
Intergrated Manufacturing, 23(11), 992-1001.