請問怎麼用matlab 畫一個傾斜的橢球？

01-08

已知橢球的三個軸長，以及每個軸的方向角，就是與X,Y,Z軸的夾角，還有橢球的中心，怎麼畫出該橢球？我知道用matlab 的函數ellipsoid(xc,yc,zc,xr,yr,zr,n)可以畫，其中：

xc,yc,zc分別表示橢球中心的x，y，z坐標。

xr,yr,zr分別表示橢球x，y，z半軸的長度。

但是該函數畫出來的橢球是標準的，請問怎麼用matlab 畫一個傾斜的橢球，就是旋轉一個標準的橢球。數據的話各位達人可以自己模擬，

設jd是3x3的方向角矩陣，其三行依次表示橢球a軸，b軸，c軸方向角，單位為弧度

則由原始坐標系旋轉到以a,b,c軸方向為坐標軸的新坐標系的旋轉矩陣顯然滿足：

M*eye(3) = cos(jd)."

顯然有

M = cos(jd)."

之後獲得變換矩陣，變換矩陣由於加上了平移變換所以是4*4的矩陣，第四行為[0 0 0 1]：

M(4,4) = 1;

第四列前三個元素為平移量，如果你需要對旋轉後的橢球進行平移（球心不在原點），那麼可以直接給M(1:3,4)賦值例如：

M(1:3,4) = 1:3; % 球心坐標為[1 2 3]

之後將旋轉矩陣輸入hgtransform即可

整合的代碼如下：

M = cos(jd)."; M(4,4) = 1; % M(1:3,4) = [xc yc zc]; % 旋轉後平移橢球心到[xc yc zc] [x,y,z] = ellipsoid(0,0,0,xr,yr,zr,n); surf(x,y,z,"Parent",hgtransform("Matrix",M)) axis image, view(3), grid on

另外想說下，這裡不建議用rotate，因為這裡通過已知的方向角求旋轉矩陣比求歐拉角方便太多（只需一步），而且rotate內部也是將求出旋轉矩陣再計算旋轉變換，所以用rotate的方法大概相當於本來你距離終點一步，結果你又跑到了起點讓rotate背著你再跑到終點。

更別說rotate在新版本中和hgtransform比起來性能先天不足的問題了

補個圖，方向角(jd)和其餘參數由隨機數生成

jd = acos(orth(rand(3))); xr = rand; yr = rand; zr = rand; n = 50;

// 冒昧來答一下這題，現學現賣，有不當之處還望大家指正。

題主問題的癥結在於如何通過一個標準的橢球(長軸和短軸分別與坐標系的 $x$ 軸和 $y$ 軸重合)經過一定的變換得到一般意義上（不那麼規則）的橢球。我在你的思路基礎上，想了一種大概可行的解決方案，但是可能偏數學一些，僅供參考。

首先，我們知道如上所說的標準橢球有以下方程唯一決定：

$frac{(x-x_0)^2}{r_x^2}+frac{(y-y_0)^2}{r_y^2}+frac{(z-z_0)^2}{r_z^2}=1$ ，其中 $(x_0, y_0, z_0)$ 表示橢球的中心， $(r_x, r_y, r_z)$ 表示各方向上對應的半軸長。

而一般意義上的橢球方程，其解析式會出現 $x, y, z$ 的交叉項，沒有這麼標準。當然你如果能根據提供的軸長，中心，方位角等信息將其轉化成關於 $x,y,z$ 的三元二次方程，然後通過等值面函數

isosurface

變相實現三維隱函數的可視化（這個思路就留給大家嘗試吧。）

下面就講講另一種「曲線救國」的思路：

一個直觀的感覺就是可以通過這種標準的橢球繞中心進行某種程度的旋轉來「矯正」其方位，使其最終和目標橢球完全重合，恰好Matalb提供繞指定軸將圖形旋轉一定角度的函數：

rotate

其使用方法還請參考官方文檔，或者移步相關問題搬運：

matlab中如何實現三維圖像以指定角度旋轉,坐標軸不變? - 知乎用戶的回答

現在問題的關鍵就是該怎麼來指定標準橢球的旋轉方式，《經典力學》中通過「歐拉角」來定義剛體的旋轉：

萊昂哈德·歐拉用歐拉角來描述剛體在三維歐幾里得空間的取向。對於任何參考系，一個剛體的取向，是依照順序，從這參考系，做三個歐拉角的旋轉而設定的。所以，剛體的取向可以用三個基本旋轉矩陣來決定。換句話說，任何關於剛體旋轉的旋轉矩陣是由三個基本旋轉矩陣複合而成的。--摘自Wikipedia

這裡可以簡單將一般意義上的橢球理解為一個剛體，其「取向」，或者說其各軸朝向可以通過「歐拉角」來唯一指定，知道其在參考系（我們平常說的空間直角坐標系，保持不變）中的歐拉角，就能通過對參考系的旋轉變換得到目標橢球的各個朝向，而標準橢球的朝向恰好又和參考系（空間直角坐標系）完全一致，意味著可以實現通過對標準橢球的旋轉來得到目標橢球。

上面闡述了「旋轉」的大致思路及可行性，下面簡單看下歐拉角的定義：

對於在三維空間里的一個參考系，任何坐標系的取向，都可以用三個歐拉角來表現。參考系又稱為實驗室參考系，是靜止不動的。而坐標系則固定於剛體，隨著剛體的旋轉而旋轉。

這裡說的參考系就是我們將的空間直角坐標系，其保持不變，而坐標系可以理解成由橢球三個軸組成的，其隨著橢球的旋轉而發生聯動。

圖片來自維基百科

三個歐拉角： $(alpha, eta, gamma)$ 。藍色的軸是 $xyz$ -軸，紅色的軸是 $XYZ$ -坐標軸。綠色的線是交點線 ( $vec{N}$ )。
參閱上圖。設定 $xyz$ -軸為參考系（空間直角坐標系）的參考軸。稱 $xy$ -平面與 $XY$ -平面的相交為交點線，用英文字母（ $vec{N}$ ）代表。
$zxz$ 順規的歐拉角可以靜態地這樣定義：

$alpha$ 是 $x$ -軸與交點線的夾角

$eta$ 是 $z$ -軸與 $Z$ -軸的夾角

$gamma$ 是交點線與 $X$ -軸的夾角

其它相關細節還請參見歐拉角。

所以我們只要能根據你提供的目標橢球的三個軸的方位角來算出這樣的歐拉角 $(alpha, eta, gamma)$ ，就能通過對標準橢球旋轉來還原出目標橢球。

為了統一表述，這裡對方位角做出如下定義：

一個向量的方位角由其與 $z$ 軸夾角 $heta$ ，其在 $xy$ -平面投影與 $x$ 軸夾角 $psi$ 唯一確定

我們假設目標橢球其坐標系 $x,z$ 軸的方位角為 $( heta_x, psi_x), ( heta_z, psi_z)$ （這裡之所以不列出 $y$ 軸方位角是因為前兩者指定後，其被唯一確定，整個計算過程也並不涉及它，所以可以不用理會）。只要一點空間幾何的知識就能算出:

$x$ 軸的方向向量 $(sin heta_xcospsi_x, sin heta_xsin psi_x, cos heta_x)$
$z$ 軸的方向向量 $(sin heta_zcospsi_z, sin heta_zsin psi_z, cos heta_z)$

然後為了計算歐拉角，還需要確定 $vec{N}$ ，考慮到其為兩個平面的交線，直接根據 $z imes Z$ 得到

$vec{N} = (-sin heta_zsin psi_z, sin heta_z cospsi_z, 0)$

OK，現在可以依次計算 $alpha, eta, gamma$ ，這裡省去計算過程，直接給出最終結果：

$alpha = frac{pi}{2} + psi_z$

$eta = heta_z$

$gamma = arccos ig(sin heta_xsin(psi_x-psi_z)ig)$

注意各角度取值範圍 $alpha, gamma in [0, 2pi], eta in [0, pi]$ 。到這裡基本已經完成了關鍵計算了，然後我們按照如下順序對標準橢球進行旋轉即可還原出目標橢球：

設 $f$ 為你根據ellipsoid計算出的標準橢球 $frac{(x-x_0)^2}{r_x^2}+frac{(y-y_0)^2}{r_y^2}+frac{(z-z_0)^2}{r_z^2}=1$ 經surface得到的圖像句柄，然後將其用rotate函數進行如下旋轉：

繞 $z$ -軸旋轉 $alpha$ 角度
繞 $x$ -軸旋轉 $eta$ 角度
繞 $z$ -軸旋轉 $gamma$ 角度

這裡採用的是 $zxz$ 旋轉，這和你按何種順序定義歐拉角有關：

在經典力學裡，時常用 $zxz$ 順規來設定歐拉角；照著第二個轉動軸的軸名，簡稱為 $x$ 順規。另外，還有別種歐拉角組。合法的歐拉角組中，唯一的限制是，任何兩個連續的旋轉，必須繞著不同的轉動軸旋轉。因此，一共有12種順規。例如， $y$ 順規，第二個轉動軸是 $y$ -軸，時常用在量子力學，核子物理學，粒子物理學。另外，還有一種順規， $xyz$ 順規，是用在航空航天工程學；參閱Tait-Bryan angles。

最後應該就能得到你要求的「斜著的」目標橢球了。但是這種方案只是分布拆解了利用歐拉角進行空間坐標系變換的主體思想，實際操作上不一定高效。

事實上，坐標系的變換本質上還是「旋轉矩陣」，如 @Falccm指出的那樣，這裡通過三次rotate來實現對標準橢球的旋轉相對很低效，實際上其相當於分步複合了三個旋轉矩陣（繞 $z$ 軸旋轉 $alpha$ ，繞 $x$ 旋轉 $eta$ ，以及最後繞 $z$ 軸旋轉 $gamma$ ）。而在事先已經知道了目標橢球三軸各方向角的情形下，完全沒有必要這樣折騰，可以直接寫出「旋轉矩陣」。

簡單來看下「旋轉矩陣」的意思（我也正好複習一下-_-||）

對於二維平面，若將坐標軸繞原點逆時針旋轉 $heta$ 角度，則新的坐標系可以這樣表示：

假設藍色線表示原始坐標系 $(e_1, e_2)$ 為單位向量，現將其逆時針旋轉 $heta$ 得到新的坐標系（紅色線） $(e_1^{$ ，不難得到：

$egin{pmatrix} e_1^{$

我們假設有一點在原始坐標系中的坐標為 $(x, y)$ ，顯然其繞原點旋轉後在新坐標系中的坐標不變，我們假設其旋轉後在原始坐標系中的坐標為 $(x^{$ ，自然就有

$(x, y) egin{pmatrix} e_1^{$ ，替換新坐標系後就能得到：

$egin{pmatrix} x^{$

我們稱 $egin{bmatrix} cos heta -sin heta\ sin heta cos heta end{bmatrix}$ 為對應的旋轉矩陣。

將其拓展到三維空間，也不難理解，繞不同的做標準旋轉對應的旋轉矩陣如下：

摘自維基百科--旋轉矩陣（點擊放大）

空間中任意兩個坐標系之間的轉換都能通過這三個基本旋轉矩陣複合實現（歐拉角的順規旋轉本質上是該方法的一個特殊情形），理解到這基本就差不多了，你只需根據目標橢球各軸的方向角求出對應的旋轉矩陣，然後對標準橢球數據進行相關旋轉就OK了，代碼參考Falccm的答案就行。

按慣例，最後舉個簡單例子，假設目標橢球的球心及各半軸長為：

[xc, yc, zc] = deal(1, 1, 1); [xr, yr, zr] = deal(6, 2, 3); [x, y, z] = ellipsoid(xc, yc, zc, xr, yr, zr, 35);

surf(x, y, z) % 顯示標準橢球狀態 axis equal

現在我們將該標準橢圓繞其中心旋轉，使得 $a, b, c$ 軸的方位角（按行排列）分別為

$egin{bmatrix} -frac{pi}{4} -frac{3}{4} pi frac{pi}{2}\ arccos(frac{1}{sqrt{3}}) arccos(frac{1}{sqrt{3}}) arccos(frac{1}{sqrt{3}})\ arccos(-frac{1}{sqrt{6}}) arccos(-frac{1}{sqrt{6}}) arccos(frac{2}{sqrt{6}}) end{bmatrix}$

「注」因不太清楚題主的方位角所謂何意，這裡指對應軸與 $x, y, z$ 軸的夾角

則相應的旋轉矩陣為：

$egin{bmatrix} frac{1}{sqrt{2}} frac{1}{sqrt{3}} frac{-1}{sqrt{6}} 1\ frac{-1}{sqrt{2}} frac{1}{sqrt{3}} frac{-1}{sqrt{6}} 1\ 0 frac{1}{sqrt{3}} frac{2}{sqrt{6}} 1\ 0 0 0 1 end{bmatrix}$ ，其相當於一個旋轉矩陣和一個平移矩陣的複合。

最後在代碼中通過指定該旋轉矩陣即可surf：

[xc, yc, zc] = deal(1, 1, 1); [xr, yr, zr] = deal(6, 2, 3); [x, y, z] = ellipsoid(xc, yc, zc, xr, yr, zr, 35);


rotate_matrix = [1/sqrt(2) -1/sqrt(2) 0;

                 1/sqrt(3) 1/sqrt(3) 1/sqrt(3);

                 -1/sqrt(6) -1/sqrt(6) 2/sqrt(6)

                 ];

rotate_matrix = [rotate_matrix.", [1 1 1]";

                0 0 0 1];

surf(x, y, z, "Parent", hgtransform("Matrix", rotate_matrix)) grid on axis equal

// the end

個人認為先根據橢球各軸長得到一個標準橢球的相應坐標，然後再將三個軸上的坐標旋轉（乘以正交矩陣）得到新的坐標然後繪圖

我來發布實現的代碼：

%% 將一個橢圓進行傾斜 rotax=20; %繞X軸旋轉的角度 rotay=40; % 繞Z軸旋轉的角度 rotaz=60; % 繞Y軸旋轉的角度 % rotamatrix 坐標系旋轉矩陣 rotamatrix=[cosd(rotaz)*cosd(rotay),sind(rotaz),-cosd(rotaz)*sind(rotay);... -sind(rotaz)*cosd(rotay)*cosd(rotax)+sind(rotay)*sind(rotax),... cosd(rotaz)*cosd(rotax),... sind(rotaz)*sind(rotay)*cosd(rotax)+cosd(rotay)*sind(rotax);... sind(rotaz)*cosd(rotay)*sind(rotax)+sind(rotay)*cosd(rotax),... -cosd(rotaz)*sind(rotax),... -sind(rotaz)*sind(rotay)*sind(rotax)+cosd(rotay)*cosd(rotax)]; %% 生成橢圓的坐標矩陣 [x,y,z] = ellipsoid(0,0,0,10,6,4,50); figure hold on subplot(2,1,1) surf(x, y, z) axis equal xlabel("X"),ylabel("Y"),zlabel("Z") % 對坐標點進行變換 coor = [x(:),y(:),z(:)]; % 網格坐標數據轉化為列坐標數據 % 數據進行矩陣變換 newcoor = coor*rotamatrix; newx = reshape(newcoor(:,1),size(x)) newy = reshape(newcoor(:,2),size(y)) newz = reshape(newcoor(:,3),size(z)) subplot(2,1,2) surf(newx, newy, newz) axis equal xlabel("X"),ylabel("Y"),zlabel("Z") hold off

結果圖如下：

view

換個視角可以嗎？

你是說旋轉標準橢球後得到的橢球么？最好給出一組數據，否則不太明白你的意思

旋轉有的很多種定義參數方式，同樣是三個轉角，所繞的軸不同，順序不同就可以了產生不同的旋轉，而且沒有哪種公認的方法，一般學術討論的時候都會事先說明採用的定義方法

「每個軸的方位角，就是與X,Y,Z軸的夾角」並不是明確的定義

你可以看下

Rotation matrix

Euler angles

先明確自己是哪種參數定義形式

而且估計把旋轉矩陣搞出來了圖你自己就會畫了

反之，如果現在稀里糊塗就隨便給個定義畫了，結果參數的定義方式和你的參數不符，畫出來的圖也不是你要的