Lập trình chuyển hệ tọa độ OCS sang WCS

Chuyển từ tọa độ OCS sang WCS từ thông tin trong DXF
Trong tệp DXF của Autocad, các đối tượng cơ bản như CIRCLE, ARC, LINE, POLYLINE, LWPOLYLINE,... ngoài tham số tọa độ (Với Group Code là 10, 20, 30) còn có thông tin khác đó là Extrusion direction (Với Group Code là 210, 220, 230). Nếu có thêm thông số extrusion direction thì chứng tỏ các giá trị tọa trên trong hệ tọa độ OCS (Object Coordinate System), bạn cần phải chuyển nó sang tọa độ thực WCS(World Coordinate System). Trong VBA, đơn giản bạn có thể dùng hàm TranslateCoordinates trong thư viện của Autocad.

Ở đây, tôi nói chi tiết hơn về thuật toán để chuyển từ tọa độ OCS sang WCS ở trong Autocad.
Đầu tiên, chúng ta xác định được hướng theo trục X và Y trong mặt phẳng XY. Autocad thực hiện tính toán này bằng thuật toán "Arbitrary Axis Algorithm" như sau:
- Vector chuẩn đầu vào, gọi là vector N
- Vector theo trục Y: Wy = (0, 1, 0)
Vector theo trục Z: Wz = (0, 0, 1)
- Bây giờ chúng ta tính vector theo trục X (arbitrary X axis) và trục Y (arbitrary Y axis) từ vector N này. Chúng được gọi là Ax và Ay, và vector N cũng được gọi là Az (arbitrary Z axis).
MÃ: CHỌN TẤT CẢ
// Tính Ax
If (abs (Nx) < 1/64) and (abs (Ny) < 1/64) then
Ax = Wy X N (Trong đó "X" là phép tích có hướng của hai vector).
Otherwise,
Ax = Wz X N.

// Chuẩn hóa Ax
Ax = Ax/|Ax| (Trong đó |Ax| là độ dài vector Ax)

// Tính Ay
Ay = N X Ax

// Chuẩn hóa Ay
Ay = Ay/|Ay|

Phép tích có hướng (cross product) hai vector được tính như sau:
MÃ: CHỌN TẤT CẢ
a = (a1, a2, a3)
b = (b1, b2, b3)
a × b = (a2b3 − a3b2, a3b1 − a1b3, a1b2 − a2b1)

Khi đó Ax, Ay, Az là vector đơn vị trong hệ tọa độ tọa OCS. Và giá trị bạn đọc trong DXF chính là giá trị trong hệ tọa độ này. Nếu coi Ax, Ay,Az là các ma trận 1 hàng 3 cột thì khi đó ma trận chuyển từ tọa độ WCS sang tọa độ OCS sẽ là:
A = (Ax, Ay, Az)

Vì Ax, Ay, Az là các vector chuẩn và đôi một vuông góc với nhau, nên ma trận chuyển từ tọa độ OCS sang WCS sẽ là:
B = A(T) = (AxT, AyT,AzT)
Bây giờ để chuyển tọa độ, đơn giản chỉ là phép nhân ma trận.

Ví dụ tính ma trận chuyển:
Giả sử ta có vector chuẩn (đọc từ DXF):
N = (-0.0539672421471262, 0.3086955872369724, 0.9496286491020867)
Tính toán ta được:
Ax = Wz X N = (-0.98505997510014909, -0.17221162985029148, 0.00000000000000000)
Ay = N X Ax = (0.16353709741440092, -0.93544117343888977, 0.31337745420585988)
Az = N = (-0.0539672421471262, 0.3086955872369724, 0.9496286491020867)
Đoạn mã nhỏ thực hiện chuyển tọa độ
Khi bạn đọc từ DXF, bạn sẽ lấy được:
- Vector chuẩn N (Group code là 210, 220, 230)
- Giá trị Elevation (Ví dụ với LWPolyline thì group code là 38). Đây chính là tọa độ z trong hệ tọa độ OCS

Trong đoạn mã dưới tôi sử dụng hai cấu trúc giống nhau đó là MYPOINT và MYVECTOR, chỉ đơn giản là để phân biệt mà thôi. HàmCalculateCrossProduct dùng để tính tích có hướng hai vector.

Code:

MÃ: CHỌN TẤT CẢ
typedef struct _MYPOINT
{
  double x;
  double y;
  double z;

  _MYPOINT()
  {
      x = 0;
      y = 0;
      z = 0;
  }

  _MYPOINT(double xx, double yy, double zz = 0)
  {
      x = xx;
      y = yy;
      z = zz;
  }
} MYPOINT;

typedef struct _MYVECTOR
{
  double x;
  double y;
  double z;
 
  _MYVECTOR()
  {
      x = 0;
      y = 0;
      z = 0;
  }

  _MYVECTOR(double xx, double yy, double zz)
  {
      x = xx;
      y = yy;
      z = zz;
  }
} MYVECTOR;

MYVECTOR CalculateCrossProduct(const MYVECTOR& v1, const MYVECTOR& v2)
{
  MYVECTOR vRet;

  // Tính tích có hướng
  vRet.x = v1.y*v2.z - v1.z*v2.y;
  vRet.y = v1.z*v2.x - v1.x*v2.z;
  vRet.z = v1.x*v2.y - v1.y*v2.x;

  // Chuẩn hóa thành vector đơn vị
  double vLen = sqrt(vRet.x*vRet.x + vRet.y*vRet.y + vRet.z*vRet.z);
  if (vLen>0)
  {
      vRet.x = vRet.x/vLen;
      vRet.y = vRet.y/vLen;
      vRet.z = vRet.z/vLen;
  }

  return vRet;
}

// Đoạn mã thực hiện chuyển
MYPOINT pOCSPts[2] = {MYPOINT(-77498.95681699999, -145452.798077, dElevation), MYPOINT(-77498.95681699999, -145456.840073, dElevation)};
MYVECTOR N = MYVECTOR(-0.0539672421471262, 0.3086955872369724, 0.9496286491020867);

MYVECTOR Wy = MYVECTOR(0, 1, 0);
MYVECTOR Wz = MYVECTOR(0, 0, 1);
MYVECTOR Ax;
MYVECTOR Ay;
if (abs(N.x)<1.0/64 && abs(N.y)<1.0/64)
  Ax = CalculateCrossProduct(Wy, N);
else
  Ax = CalculateCrossProduct(Wz, N);
Ay = CalculateCrossProduct(N, Ax);
MYVECTOR Az = N;

MYPOINT pWCSPts[2];
for (int i=0; i<2; i++)
{
  pWCSPts[i].x = Ax.x*pOCSPts[i].x + Ay.x*pOCSPts[i].y + Az.x*pOCSPts[i].z;
  pWCSPts[i].y = Ax.y*pOCSPts[i].x + Ay.y*pOCSPts[i].y + Az.y*pOCSPts[i].z;
  pWCSPts[i].z = Ax.z*pOCSPts[i].x + Ay.z*pOCSPts[i].y + Az.z*pOCSPts[i].z;
}


pOCSPts là tọa độ OCS và pWCSPts là tọa độ WCS. Sau khi tính nếu tọa độ z bạn không dùng bạn có thể bỏ qua.

Ở đây tôi cũng nói thêm việc chuyển hệ tọa độ này cũng được thực hiện bởi ngôn ngữ AutoLISP các bạn có thể tham khảo cuốn "hướng dẫn lập trình AutoLISP và Visual LISP tập I" của Nguyễn Hữu Lộc và Nguyễn Thanh Trung(2003). Nếu ai ở Đại Học Nha Trang có thể liên hệ trực tiếp với tôi, tôi sẽ cho mượn cuốn này