在两个任意定向的椭球之间进行插值 [英] Interpolating between two arbitrarily oriented ellipsoids

问题描述

我不确定在 Mathoverflow 上问这个问题是否更好，但我想我会先在这里查看.我尽量做到简洁明了；如果有什么需要澄清的，请告诉我.

背景

我在 R3 中有两组点，它们以(或多或少)任意方向的椭圆体的形式分布.我希望在这两个椭球之间插入一个管状结构.我也有这个管状结构所需中心线的坐标.

我使用在 Matlab 中实现的 Khachiyan 算法[1] 使用包含椭圆体的最小体积来近似椭圆体的两端，该算法返回椭圆体中心的坐标 (C) 和中心形式的椭圆矩阵 (A)，使得:

(x - C)' * A * (x - C) = 1

然后我使用奇异值分解提取椭圆体的轴长度 (a,b,c) 和旋转矩阵 (V):

[U,D,V] = svd(A);a = 1/sqrt(D(1,1));b = 1/sqrt(D(2,2));c = 1/sqrt(D(3,3));

我可以轻松插入轴长度参数(例如线性、样条).为了在方向之间进行插值，我首先将旋转矩阵转换为四元数表示.然后对于沿中心线的每个点，我使用在另一个 Matlab 文件中实现的球面线性插值 (SLERP) [2]:

for iPoint = 1 : nPointst = iPoint/(nPoints + 2);quat = slerp(startQuat,endQuat,t,0.001);R = quat2rot(quat);结尾

这就是我卡住的地方.

不幸的是，即使 SLERP在其四元数端点之间给出了一条最直和最短的路径"，[3] 得到的内插椭球有时在错误"方向旋转.也就是说，内插不是产生光滑的管子，而是产生一种扭曲的椭圆圆柱体(见下图).

我尝试检查两个四元数的点积是否为负，如果是，则使用 quatinv 反转其中一个.但是，反转会导致完全错误的结果(请参见下面的第二张附图).

我的问题是:为什么会发生这种情况，我可以做些什么来纠正这种行为?也就是说，如何沿两个椭球方向之间的真实"最短路径进行插值?

任何建议将不胜感激！

更新

我已经创建了一个最小的工作示例和一个必需的数据文件.我还附上了结果的截图.我已经将它们压缩并上传到 Dropbox.[4]

[1] 发布了一个新问题.

I'm not sure if this would be better asked on Mathoverflow, but I thought I would check here first. I have tried to be as clear and concise as possible; if there is anything that needs clearing up please let me know.

Background

I have two sets of points in R3 that are distributed in the form of (more-or-less) arbitrarily oriented ellipsoids. I wish to interpolate a tubular structure between these two ellipsoids. I also have coordinates of the desired centre line of this tubular structure.

I approximate the ellipsoids at either end with a minimum volume enclosing ellipsoid using the Khachiyan Algorithm implemented in Matlab, [1] which returns the coordinates of the centre of the ellipsoid (C), and matrix of the ellipse in centre form (A), such that:

(x - C)' * A * (x - C) = 1 

I then extract the ellipsoid's axes lengths (a,b,c) and the rotation matrix (V) using singular value decomposition:

[U,D,V] = svd(A);
a = 1/sqrt(D(1,1));
b = 1/sqrt(D(2,2));
c = 1/sqrt(D(3,3));

I can easily interpolate the axes length parameters (e.g. linear, spline). To interpolate between the orientations, I first convert the rotation matrices to quaternion representation. Then for each point along the centre line, I use spherical linear interpolation (SLERP) implemented in another Matlab file [2]:

for iPoint = 1 : nPoints
  t = iPoint / (nPoints + 2);
  quat = slerp(startQuat,endQuat,t,0.001);
  R = quat2rot(quat);
end

This is where I get stuck.

Unfortunately, even though SLERP "gives a straightest and shortest path between its quaternion endpoints," [3] the resulting interpolated ellipsoids are sometimes rotating in the "wrong" direction. That is, rather than resulting in a smooth tube, the interpolation results in a sort of twisted elliptical cylinder (see attached image, below).

I have tried checking to see if the dot product of the two quaternions is negative and if so, inverting one of them using quatinv. However, inverting results in something completely incorrect (see second attached image, below).

My question is: why is this happening, and what can I do to correct for this behavior? That is, how can I interpolate along the "true" shortest path between the two ellipsoid orientations?

Any suggestions would be greatly appreciated!

UPDATE

I have created a minimum working example and a required data file. I have also attached a screenshot of the result. I've zipped these up and uploaded them to Dropbox. [4]

[1] http://www.mathworks.com/matlabcentral/fileexchange/9542-minimum-volume-enclosing-ellipsoid/content/MinVolEllipse.m

[2] http://www.mathworks.com/matlabcentral/fileexchange/11827-slerp/content/slerp.m

[3] https://en.wikipedia.org/wiki/Slerp

[4] https://dl.dropboxusercontent.com/u/38218/ellipsoidInterpolation.zip

解决方案

The solution was to rotate everything by the inverse of the rotation matrix of one of the reference ellipsoids, such that that reference ellipsoid was axis-aligned (i.e. had no rotation). Then after interpolating each ellipsoid, rotate it back to the original reference frame by multiplying by the original rotation matrix.

I've attached a screenshot of the result:

Update

Apparently this does not work in every case. I have posted a new question here.

这篇关于在两个任意定向的椭球之间进行插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！