球坐标中的插值曲线
Interpolating curve in spherical coordinates
我有点,属于3D中的一些非参数曲线space。
我需要在球坐标系中工作,并在给定的 azimuth
和 elevation
值中插入给定的曲线。
你能帮我看看这种插值算法吗?
P.S. 我不确定,从我的问题想法 1D 曲线绘制在 1D space 很清楚。因此,以防万一下面的代码创建了这样的参数。
x=linspace(-2*pi,2*pi,10^3);
y=sin(x);
z=sinh(x);
更新:
感谢@jodag 非常好的 post 关于上采样。 slerp
完美!尽管如此,我的问题是不同的,或者我不明白如何使用上采样来解决它。
假设我将我的点定义为数组(只是一些用于演示的虚拟数据):
x_given = rand(9, 1) - 0.5;
y_given = rand(9, 1) - 0.5;
z_given = rand(9, 1) - 0.5;
x_given(end+1) = x_given(1); % should work for closed curves
y_given(end+1) = y_given(1); % should work for closed curves
z_given(end+1) = x_given(1); % should work for closed curves
我想找到给定 elevation_q
和 azimuth_q
的所有点 (x_i, y_i, z_i
):
elevation_q = pi./4;
azimuth_q = pi./7;
[x_i, y_i ,z_i] = interpolateCurve(x_given, y_given, y_given, elevation_q, azimuth_q);
那是我真正需要的。我明白,我们不需要从笛卡尔坐标系来回跳到球坐标系。但即使我将 x_given, y_given, y_given
替换为 elevation_given, azimuth_given, r_given
,我也不知道如何编写 interpolateCurve()
函数。
这是一个在更高维度中对曲线进行线性插值的示例 space。请注意,我引入了一个参数化变量 t
对应于点的索引。
function [xq,yq,zq] = interp_sph(x,y,z,upsample)
[az,el,r] = cart2sph(x,y,z);
t = 1:numel(x);
tq = 1:(1/upsample):numel(x);
azq = interp1(t,az,tq);
elq = interp1(t,el,tq);
rq = interp1(t,r,tq);
[xq,yq,zq] = sph2cart(azq,elq,rq);
end
从问题中不清楚这是否是您真正要查找的内容。球坐标中的直接插值肯定存在问题。例如,如果两个点彼此靠近,但一个点的方位角 = -179 度,另一个点的方位角 = 179,则插值将产生跨越 [-179,179] 的样本,而不是采用最短的 2 度路径。您可能感兴趣的是 slerp,它给出了更好的插值结果。
我写了下面的slerp例子来演示。请注意,slerp 要求定义 1/sin(theta)
,其中 theta
是函数中两点之间的角度。如果 theta == pi
那么我们有问题,slerp 将无法工作,所以如果您的数据接近原点,请小心。
function [xq,yq,zq] = interp_slerp(x,y,z,upsample)
r = sqrt(x.^2+y.^2+z.^2);
vq = zeros(3,(upsample-1)*(numel(x)-1)+1);
for idx = 1:numel(x)-1
t = 0:1;
tq = linspace(0,1,upsample);
% compute interpolated direction
dq = slerp([x(idx),y(idx),z(idx)],[x(idx+1),y(idx+1),z(idx+1)],tq);
% linearly interpolate radius
rq = interp1(t,r(idx:(idx+1)),tq);
% get interpolated path
idxq = (idx-1)*upsample + 1;
vq(:,idxq:(idxq+(upsample-1))) = bsxfun(@times,rq,dq);
end
xq = vq(1,:);
yq = vq(2,:);
zq = vq(3,:);
end
function pq = slerp(p0,p1,t)
p0 = p0(:) / norm(p0);
p1 = p1(:) / norm(p1);
theta = acos( dot(p0,p1) );
if(0==theta)
pq = a;
elseif(pi==theta)
error('Angle between points cannot be exactly pi.');
else
pq = bsxfun(@times,(sin((1.0-t)*theta)/sin(theta)),p0) + bsxfun(@times,(sin(t*theta)/sin(theta)),p1);
end
end
为了演示 interp_sph
的问题,请考虑以下示例。
% create test path
az = pi+sin(linspace(pi/2,3*pi/2,10));
el = cos(linspace(-pi,pi,10));
r = 1 + 1/5*sin(6*az).*sin(5*el);
[x,y,z] = sph2cart(az,el,r);
% test both spherical interpolation and slerp
upsample = 20;
[xq_sph,yq_sph,zq_sph] = interp_sph(x,y,z,upsample);
[xq_slerp,yq_slerp,zq_slerp] = interp_slerp(x,y,z,upsample);
% plot results
plot3(x,y,z,'LineWidth',2); hold on;
plot3(xq_sph,yq_sph,zq_sph,'LineWidth',2);
plot3(xq_slerp,yq_slerp,zq_slerp,'LineWidth',2);
axis('vis3d');
axis([-1,1,-1,1,-1,1]);
grid('on');
legend('Cartisian', 'Spherical', 'SLERP');
这产生了下面的情节。请注意,当方位角从正变为负时,球面插值会走很长一段路,但 slerp 会走预期路线。
我有点,属于3D中的一些非参数曲线space。
我需要在球坐标系中工作,并在给定的 azimuth
和 elevation
值中插入给定的曲线。
你能帮我看看这种插值算法吗?
P.S. 我不确定,从我的问题想法 1D 曲线绘制在 1D space 很清楚。因此,以防万一下面的代码创建了这样的参数。
x=linspace(-2*pi,2*pi,10^3);
y=sin(x);
z=sinh(x);
更新:
感谢@jodag 非常好的 post 关于上采样。 slerp
完美!尽管如此,我的问题是不同的,或者我不明白如何使用上采样来解决它。
假设我将我的点定义为数组(只是一些用于演示的虚拟数据):
x_given = rand(9, 1) - 0.5;
y_given = rand(9, 1) - 0.5;
z_given = rand(9, 1) - 0.5;
x_given(end+1) = x_given(1); % should work for closed curves
y_given(end+1) = y_given(1); % should work for closed curves
z_given(end+1) = x_given(1); % should work for closed curves
我想找到给定 elevation_q
和 azimuth_q
的所有点 (x_i, y_i, z_i
):
elevation_q = pi./4;
azimuth_q = pi./7;
[x_i, y_i ,z_i] = interpolateCurve(x_given, y_given, y_given, elevation_q, azimuth_q);
那是我真正需要的。我明白,我们不需要从笛卡尔坐标系来回跳到球坐标系。但即使我将 x_given, y_given, y_given
替换为 elevation_given, azimuth_given, r_given
,我也不知道如何编写 interpolateCurve()
函数。
这是一个在更高维度中对曲线进行线性插值的示例 space。请注意,我引入了一个参数化变量 t
对应于点的索引。
function [xq,yq,zq] = interp_sph(x,y,z,upsample)
[az,el,r] = cart2sph(x,y,z);
t = 1:numel(x);
tq = 1:(1/upsample):numel(x);
azq = interp1(t,az,tq);
elq = interp1(t,el,tq);
rq = interp1(t,r,tq);
[xq,yq,zq] = sph2cart(azq,elq,rq);
end
从问题中不清楚这是否是您真正要查找的内容。球坐标中的直接插值肯定存在问题。例如,如果两个点彼此靠近,但一个点的方位角 = -179 度,另一个点的方位角 = 179,则插值将产生跨越 [-179,179] 的样本,而不是采用最短的 2 度路径。您可能感兴趣的是 slerp,它给出了更好的插值结果。
我写了下面的slerp例子来演示。请注意,slerp 要求定义 1/sin(theta)
,其中 theta
是函数中两点之间的角度。如果 theta == pi
那么我们有问题,slerp 将无法工作,所以如果您的数据接近原点,请小心。
function [xq,yq,zq] = interp_slerp(x,y,z,upsample)
r = sqrt(x.^2+y.^2+z.^2);
vq = zeros(3,(upsample-1)*(numel(x)-1)+1);
for idx = 1:numel(x)-1
t = 0:1;
tq = linspace(0,1,upsample);
% compute interpolated direction
dq = slerp([x(idx),y(idx),z(idx)],[x(idx+1),y(idx+1),z(idx+1)],tq);
% linearly interpolate radius
rq = interp1(t,r(idx:(idx+1)),tq);
% get interpolated path
idxq = (idx-1)*upsample + 1;
vq(:,idxq:(idxq+(upsample-1))) = bsxfun(@times,rq,dq);
end
xq = vq(1,:);
yq = vq(2,:);
zq = vq(3,:);
end
function pq = slerp(p0,p1,t)
p0 = p0(:) / norm(p0);
p1 = p1(:) / norm(p1);
theta = acos( dot(p0,p1) );
if(0==theta)
pq = a;
elseif(pi==theta)
error('Angle between points cannot be exactly pi.');
else
pq = bsxfun(@times,(sin((1.0-t)*theta)/sin(theta)),p0) + bsxfun(@times,(sin(t*theta)/sin(theta)),p1);
end
end
为了演示 interp_sph
的问题,请考虑以下示例。
% create test path
az = pi+sin(linspace(pi/2,3*pi/2,10));
el = cos(linspace(-pi,pi,10));
r = 1 + 1/5*sin(6*az).*sin(5*el);
[x,y,z] = sph2cart(az,el,r);
% test both spherical interpolation and slerp
upsample = 20;
[xq_sph,yq_sph,zq_sph] = interp_sph(x,y,z,upsample);
[xq_slerp,yq_slerp,zq_slerp] = interp_slerp(x,y,z,upsample);
% plot results
plot3(x,y,z,'LineWidth',2); hold on;
plot3(xq_sph,yq_sph,zq_sph,'LineWidth',2);
plot3(xq_slerp,yq_slerp,zq_slerp,'LineWidth',2);
axis('vis3d');
axis([-1,1,-1,1,-1,1]);
grid('on');
legend('Cartisian', 'Spherical', 'SLERP');
这产生了下面的情节。请注意,当方位角从正变为负时,球面插值会走很长一段路,但 slerp 会走预期路线。