在不同大小的球体上贴片

Patch on a sphere of varying size

想象一下粘在球体上的补丁。当 球体按比例放大或缩小 时,我如何设法使补丁保持其中心位置和表面积?通常,只有补丁的曲率应该改变,因为它是“粘”到球体上的。假设补丁被描述为一组(纬度,经度)坐标。

一种可能的解决方案包括将补丁的地理坐标转换为地心坐标(从上方直接垂直观察的补丁),从而制作 2D 纹理,然后随着球体大小的变化而放大或缩小。但我不确定这是否是正确的方法,以及这与预期效果的接近程度。

我是新手,所以也许 Unity 可以在应用纹理时使用正确的设置选项简单地做到这一点。在这种情况下,应该为纹理使用哪个输入贴图投影?或者也许我应该使用 3D 表面并以某种方式将其“钉”到球体上。

谢谢!!

编辑

我正在添加一个插图来展示当球体按比例放大或缩小时补丁应该如何变形。在一个非常小的球体上,补丁最终会环绕。而在更大的球体上,补丁几乎是平坦的。贴片的变形可以被认为类似于将同一张贴纸粘在不同大小的球体上。

补丁的几何形状可以是任何多边形表面,并且如前所述,在球体按比例放大或缩小时必须保持其中心位置和表面积。

你可以,例如做类似

的事情
public class Example : MonoBehaviour
{
    public Transform sphere;

    public float latitude;
    public float longitude;

    private void Update()
    {
        transform.position = sphere.position
                             + Quaternion.AngleAxis(longitude, -Vector3.up) 
                             * Quaternion.AngleAxis(latitude, -Vector3.right) 
                             * sphere.forward * sphere.lossyScale.x / 2f;
        transform.LookAt(sphere);
        transform.Rotate(90,0,0);
    }
}

图钉 不是 球体的子项。它会产生一个像这样的图钉(红色):


或者如前所述,您可以使图钉成为球体的子结构,例如

Sphere
|--PinAnchor
   |--Pin

因此,为了更改 Pin 位置,您需要旋转 PinAnchor。 Pin 本身会更新自己的比例,因此它始终具有特定的目标比例,例如喜欢

public class Example : MonoBehaviour
{
    public float targetScale;

    private void Update()
    {
        var scale = transform.parent.lossyScale;
        var invertScale = new Vector3(1 / scale.x, 1 / scale.y, 1 / scale.z);
        if (float.IsNaN(invertScale.x)) invertScale.x = 0;
        if (float.IsNaN(invertScale.y)) invertScale.y = 0;
        if (float.IsNaN(invertScale.z)) invertScale.z = 0;

        transform.localScale = invertScale * targetScale;
    }
}

假设您有一个半径为 R1 的球体,其中心位于标准坐标系 O e1 e2 e3 的原点。然后球体由 3D 中满足方程 x[0]^2 + x[1]^2 + x[2]^2 = R1^2 的所有点 x = [x[0], x[1], x[2]] 给出。在这个球体上,你有一个补丁,补丁有一个中心 c = [c[0], c[1], c[2]]。 首先,旋转面片,使中心 c 指向北极,然后将其投影到平面上,使用半径为 R1 的球体的区域保留地图,然后使用类似的方法将其映射回来区域保留地图,但半径为 R2 球体,最后将北极旋转回中心的缩放位置。

您可能需要定义的函数:

函数1:定义球面坐标

x = sc(u, v, R):
    return
        x[0] = R*sin(u)*sin(v)
        x[1] = R*sin(u)*cos(v)
        x[2] = R*cos(u)

where
0 <= u <= pi and 0 <= v < 2*pi

函数2:定义反球坐标:

[u, v] = inv_sc(x, R):
    return
        u = arccos( x[2] / R )
        if x[1] > 0
           v = arccot(x[0] / x[1]) if x[1] > 0 
        else if x[1] < 0 
           v = 2*pi - arccot(x[1] / x[0]) 
        else if x[1] = 0 and x[0] > 0
           v = 0
        else if x[1] = 0 and x[0] < 0
           v = pi

where  x[0]^2 + x[1]^2 + x[2]^2 = R^2

函数3:将中心c旋转到北极的旋转矩阵:

假设中心 c 在球坐标 [uc, vc] 中给出。然后应用函数 1

c = [c[0], c[2], c[3]] = sc(uc, vc, R1)

然后,找到我们有 c[i] = min( abs(c[0]), abs(c[1]), abs(c[2])) 的索引 i。说 i=2 并取坐标向量 e2 = [0, 1, 0].

计算叉积向量cross(c, e2)cross(cross(c, e2), c),将它们视为行向量,并形成3 by 3旋转矩阵

        A3 = c / norm(c)
        A2 = cross(c, e2) / norm(cross(c, e2))
        A1 = cross(A2, A3)
        A = [ A1,
              A2,
              A3 ]

函数 4:

[w,z] = area_pres(u,v,R1,R2):
    return
       w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )
       z = v

现在,如果您将球体从半径 R1 重新缩放到半径 R2,那么球体上半径为 R1 的补丁中的任何点 x 都会得到 t 运行通过以下 t运行 变形链形成半径 R2 球体上的点 y

If x is given in spherical coordinates `[ux, vx]`, first apply

x = [x[0], x[1], x[2]] = sc(ux, vx, R1)

Then rotate with the matrix A:

x = matrix_times_vector(A, x)

Then apply the chain of transformations:

[u,v] =  inv_sc(x, R1)
[w,z] = area_pres(u,v,R1,R2)
y = sc(w,z,R2)

Now y is on the R2 sphere. 
Finally,
 
y = matrix_times_vector(transpose(A), y)

因此所有这些点 y 填充了半径为 R2 的球体上相应的 t运行sformed patch 并且 R2 上的 patch-area 等于 patch -球体上原始补丁的区域 R1。再加上中心点 c 只是沿着从球体中心发出的光线按比例放大或缩小。

这个方法背后的总体思路是,基本上,R1 球体的面积元素是 R1^2*sin(u) du dv,我们可以寻找经纬度坐标 [=60] 的 t运行 形式=] 到 R2 球体的经纬度坐标 [w,z] 中,我们有函数 w = w(u,v)z = z(u,v) 这样

R2^2*sin(w) dw dz = R1^2*sin(u) du dv

当你展开 [w,z] 关于 [u,v] 的导数时,你得到

dw = dw/du(u,v) du + dw/dv(u,v) dv
dz = dz/du(u,v) du + dz/dv(u,v) dv

将它们代入第一个公式,你得到

R2^2*sin(w) dw dz = R2^2*sin(w) * ( dw/du(u,v) du + dw/dv(u,v) dv ) wedge ( dz/du(u,v) du + dz/dv(u,v) dv ) 
      = R1^2*sin(u) du dv

简化为等式

R2^2*sin(w) * ( dw/du(u,v) dz/dv(u,v)  -  dw/dv(u,v) dz/du(u,v) ) du dv = R^2*sin(u) du dv  

所以gua运行求出R1上的球面贴片与其在R2上的像的t运行变形的面积属性的一般微分方程为

R2^2*sin(w) * ( dw/du(u,v) dz/dv(u,v)  -  dw/dv(u,v) dz/du(u,v) ) = R^2*sin(u)

现在回想一下,补丁的中心已经旋转到 R1 球体的北极,所以您可以认为补丁的中心是北极。如果你想要补丁的 t运行 变形,使其从补丁的中心有点均匀和各向同性,即当站在补丁的中心 c 时(c = north pole)你看到变形的补丁,以便保留经度(通过 c 的大圆)(即来自经度的所有点都映射到相同经度的点),你得到经度坐标 v 的限制点 [u, v] 的 t运行 形成一个新点 [w, z],它应该在同一经度上,即 z = v。因此这样的经度保留 t运行sformation 应该是这样的:

w = w(u,v)
z = v

因此,面积保持方程简化为以下偏微分方程

R2^2*sin(w) * dw/du(u,v) = R1^2*sin(u)

因为 dz/dv = 1dz/du = 0.
求解,先固定变量v,得到常微分方程

R2^2*sin(w) * dw = R1^2*sin(u) du

谁的解

R2^2*(1 - cos(w)) = R1^2*(1 - cos(u)) + const

因此,当你让v变化时,偏微分方程的通解

R2^2*sin(w) * dw/du(u,v) = R^2*sin(u)

隐式形式(链接变量 w, u, v 的等式)应该类似于

R2^2*(1 - cos(w)) = R1^2*(1 - cos(u)) + f(v)

for any function f(v)

但是,我们不要忘记北极在这个 t运行 变化期间保持不变,即我们有限制,即 w= 0 每当 u = 0。将此条件代入上面的等式,您将得到函数 f(v)

的限制条件
R2^2*(1 - cos(0)) = R1^2*(1 - cos(0)) + f(v)
R2^2*(1 - 1) = R1^2*(1 - 1) + f(v)
0 = f(v)

for every longitude v

因此,一旦您将经度运行变形为相同的经度并保留北极,剩下的唯一选择就是等式

R2^2*(1 - cos(w)) = R1^2*(1 - cos(u))

这意味着当你解决 w 你得到

w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )

因此,球体R1上的小块与球体R2上的小块面积相同,中心固定,中心均匀变形,使得经度为t的对应面积保持t运行变形运行变形为相同的经度,是

w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )
z = v

这里我在Python和运行中实现了一些简单的模拟功能:

import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def trig(uv):
    return np.cos(uv), np.sin(uv)

def sc_trig(cos_uv, sin_uv, R):
    n, dim = cos_uv.shape
    x = np.empty((n,3), dtype=float)
    x[:,0] = sin_uv[:,0]*cos_uv[:,1] #cos_u*sin_v
    x[:,1] = sin_uv[:,0]*sin_uv[:,1] #cos_u*cos_v
    x[:,2] = cos_uv[:,0] #sin_u
    return R*x

def sc(uv,R):
    cos_uv, sin_uv = trig(uv)
    return sc_trig(cos_uv, sin_uv, R)

def inv_sc_trig(x):
    n, dim = x.shape
    cos_uv = np.empty((n,2), dtype=float)
    sin_uv = np.empty((n,2), dtype=float)
    Rad = np.sqrt(x[:,0]**2 + x[:,1]**2 + x[:,2]**2)
    r_xy = np.sqrt(x[:,0]**2 + x[:,1]**2)
    cos_uv[:,0] = x[:,2]/Rad #cos_u = x[:,2]/R
    sin_uv[:,0] = r_xy/Rad #sin_v = x[:,1]/R
    cos_uv[:,1] = x[:,0]/r_xy
    sin_uv[:,1] = x[:,1]/r_xy
    return  cos_uv, sin_uv

def center_x(x,R):
    n, dim = x.shape
    c = np.sum(x, axis=0)/n
    return R*c/math.sqrt(c.dot(c))

def center_uv(uv,R):
    x = sc(uv,R)
    return center_x(x,R)
    
def center_trig(cos_uv, sin_uv, R):
    x = sc_trig(cos_uv, sin_uv, R)
    return center_x(x,R)

def rot_mtrx(c):
    i = np.where(c == min(c))[0][0]
    e_i = np.zeros(3)
    e_i[i] = 1
    A = np.empty((3,3), dtype=float)
    A[2,:] = c/math.sqrt(c.dot(c))
    A[1,:] = np.cross(A[2,:], e_i)
    A[1,:] = A[1,:]/math.sqrt(A[1,:].dot(A[1,:]))
    A[0,:] = np.cross(A[1,:], A[2,:])
    return A.T # ready to apply to a n x 2 matrix of points from the right

def area_pres(cos_uv, sin_uv, R1, R2):
    cos_wz = np.empty(cos_uv.shape, dtype=float)
    sin_wz = np.empty(sin_uv.shape, dtype=float)
    cos_wz[:,0] = 1 - (R1/R2)**2 * (1 - cos_uv[:,0])
    cos_wz[:,1] = cos_uv[:,1]
    sin_wz[:,0] = np.sqrt(1 - cos_wz[:,0]**2)
    sin_wz[:,1] = sin_uv[:,1]
    return cos_wz, sin_wz

def sym_patch_0(n,m):    
    u = math.pi/2 + np.linspace(-math.pi/3, math.pi/3, num=n)
    v = math.pi/2 + np.linspace(-math.pi/3, math.pi/3, num=m)
    uv = np.empty((n, m, 2), dtype=float)
    uv[:,:,0] = u[:, np.newaxis]
    uv[:,:,1] = v[np.newaxis,:]
    uv = np.reshape(uv, (n*m, 2), order='F')
    return uv, u, v

uv, u, v = sym_patch_0(18,18)
r1 = 1
r2 = 2/3
r3 = 2
limits = max(r1,r2,r3)

p = math.pi

x = sc(uv,r1) 

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(x[:,0], x[:,1], x[:,2])

ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

B = rot_mtrx(center_x(x,r1))
x = x.dot(B)
cs, sn = inv_sc_trig(x)

cs1, sn1 = area_pres(cs, sn, r1, r2)
y = sc_trig(cs1, sn1, r2)
y = y.dot(B.T)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(y[:,0], y[:,1], y[:,2])

ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

cs1, sn1 = area_pres(cs, sn, r1, r3)
y = sc_trig(cs1, sn1, r3)
y = y.dot(B.T)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(y[:,0], y[:,1], y[:,2])

ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

我们可以看到当球体的半径从半径 2/3、通过半径 1 最后变为半径 2 时补丁如何变形的三个图形。补丁的面积没有改变,并且 t运行贴片的变形在所有方向上都是均匀的,没有过度变形。

我要添加另一个答案,因为您可能会认为不同的属性对您的补丁转换很重要,更具体地说,具有最小(某种意义上)的失真,并且补丁的区域保留不是一样重要。

假设您要创建从半径为 R1 的球体上的面片(具有相对良好边界的球体的开放子集,例如分段平滑或什至分段测地线边界)到半径为 R2 的球体上相应的补丁。但是,您希望转换不扭曲 R1 上的原始补丁并将其映射到 R2。假设 R1 上的补丁有一个区别点 c,称为中心。这可能是它的几何中心,即它的质心(重心),或者以其他方式选择的点。

对于此讨论,我们假设中心 c 位于球体 R1 的北极。如果不是,我们可以简单地将它旋转到北极(旋转中心的一种方法见我之前的post),这样标准的球坐标[u, v](纬度和经度)自然适用, 即

for sphere R1:
x[0] = R1*sin(u)*cos(v)
x[1] = R1*sin(u)*sin(v)
x[2] = R1*cos(u)

for sphere R2:
y[0] = R2*sin(w)*cos(z)
y[1] = R2*sin(w)*sin(z)
y[2] = R2*cos(w)

c 的坐标为 [0,0](或任何 [0,v],因为这些坐标在极点处具有奇点)。理想情况下,如果您在两个面片之间构建等距变换(等距是一种保留距离、角度和面积的变换),那么您就完成了。然而,这两个球体具有不同的半径 R1R2,因此它们具有不同的固有曲率,因此贴片之间不可能存在等距。尽管如此,让我们看看等距会做什么:等距是一种变换,它将第一个球体的度量张量(线元素,我们在球体上测量距离的方式)转换为第二个球体的度量张量,即

Metric tensor of R1:
R1^2 * ( du^2 + (sin(u))^2 dv^2 )

Metric tensor of R2: 
R2^2 * ( dw^2 + (sin(w))^2 dz^2 )

An isometry: [u,v] --> [w,z] so that
R1^2 * ( du^2 + (sin(u))^2 dv^2 ) = R2^2 * ( dw^2 + (sin(w))^2 dz^2 )

等距会做什么,首先它将球面测地线(大圆)发送到球面测地线,所以特别是 R1 的纵向圆应该映射到 R2 的纵向圆,因为我们希望 R1 的北极映射到 R2 的北极。此外,等距会保留角度,因此特别是,它会保留纵向圆之间的角度。由于零经度圆与经度经度圆 v 之间的夹角等于 v (如果加上球体绕北极的全局旋转,则达到一个常数的平移,但是我们不想要那个),那么 v 应该由等距图保存(即等轴测图应该保留北极的方位角)。这意味着补丁之间所需的等轴测图应具有

的形式
Map between patch on R1 and patch on R2, 
which maps the north pole of R1 to the north pole of R2:

w = w(u, v)
z = v

此外,由于球体在任何点和任何方向看起来都相同(它在任何地方都是均匀且各向同性的),特别是对于北极来说是这样,因此等轴测图在观察时应该在所有方向上相同地变换从北极(术语是“等距变换应该与表面的等距自同构组交换”)这产生 w = w(u, v) 不应该依赖于变量 v:

Map between patch on R1 and patch on R2, 
which maps the north pole of R1 to the north pole of R2:

w = w(u)
z = v

找到 R1R2 上的补丁之间的等距变换的最后一步是确保变换前后的度量张量相等,即:

R2^2 * ( dw^2 + (sin(w))^2 dz^2 ) = R1^2 * ( du^2 + (sin(u))^2 dv^2 )

dw = (dw/du(u)) du   and  dz = dv

R2^2 * ( (dw/du(u))^2 du^2 + (sin( w(u) ))^2 dv^2 ) = R1^2 * ( du^2 + (sin(u))^2 dv^2 )

set K = R1/R2

( dw/du(u) )^2 du^2 + (sin( w(u) ))^2 dv^2  = K^2 du^2 + K^2*(sin(u))^2 dv^2 

为了使后一个方程成立,我们需要函数w = w(u)满足以下两个限制

dw/du(u) = K

sin(w(u)) = K * sin(u)

然而,我们只有一个函数w(u)和两个只有在K = 1(即R1 = R2)时才满足的方程,事实并非如此。这就是等距条件中断的地方,这就是为什么在 R1 != R2 时球体 R1 上的补丁和 R2 上的补丁之间没有等距变换的原因。我们可以尝试做的一件事是找到一种变换,在某种合理的意义上最小化度量张量之间的差异(即我们希望以某种方式最小化变换的非等距度 [w = w(u), z = v] )。为此,我们可以定义一个拉格朗日差异函数(是的,就像在物理学中一样)并尝试将其最小化:

Lagrangian:
L(u, w, dw/du) = ( dw/du - K )^2 + ( sin(w) - K*sin(u) )^2

minimize the action: 
S[w] = integral_0^u2  L(u, w(u), dw/du(u))du  

or more explicitly, find the function `w(u)` that makes 
the sum (integral) of all discrepancies:
S[w] = integral_0^u2 ( ( dw/du(u) - K )^2 + ( sin(w(u)) - K*sin(u) )^2 )du
minimal

为了找到最小化上述差异积分S[w]的函数w(u),需要推导与拉格朗日方程L(u, w, dw,du)相关的欧拉-拉格朗日方程并求解它们.本例中的欧拉-拉格朗日方程为一且为二阶导数一:

d^2w/du^2 = sin(w)*cos(w) - K*sin(u)*cos(w)
w(0) = 0
dw/du(0) = K 

或使用替代符号:

w''(u) = sin(w(u))*cos(w(u)) - K*sin(u)*cos(w(u))
w(0) = 0
w'(0) = K

条件w'(0) = K的原因来自强加等距身份

( dw/du(u) )^2 du^2 + (sin( w(u) ))^2 dv^2  = K^2 du^2 + K^2*(sin(u))^2 dv^2 

u = 0时,我们已经知道w(0) = 0,因为我们要将北极映射到北极,所以后面的恒等式简化为

( dw/du(0) )^2 du^2 + (sin(0))^2 dv^2  = K^2 du^2 + K^2*(sin(0))^2 dv^2

( dw/du(0) )^2 du^2 = K^2 du^2 

( dw/du(0) )^2 = K^2 

时成立
dw/du(0) = u'(0) = K  

现在,为了分别在两个半径为 R1R2 的球体上获得一个北极点变换,失真尽可能小(关于拉格纳吉误差), 我们要解决非线性初值问题

d^2w/du^2 = sin(w)*cos(w) - K*sin(u)*cos(w)
w(0) = 0
dw/du(0) = K 

或写成两个一阶导数微分方程组(哈密尔顿形式):

dw/du = p
dp/du = sin(w)*cos(w) - K*sin(u)*cos(w)
w(0) = 0
p(0) = K 

我严重怀疑这是一个完全可解(可积)的常微分方程组,但是具有相当小的积分步长的数值积分可以给出一个很好的离散解,它结合了一个很好的插值方案,比如三次样条曲线,可以给你一个非常准确的解决方案。

现在,如果您不太关心补丁之间完全相等的区域,而是相当接近的区域,并且实际上更喜欢尽可能小的(在某种意义上)几何变形,您可以简单地使用这个模型和停在这里。但是,如果你真的坚持两个补丁之间的面积相等,你可以继续,通过将球体 R1 上的原始补丁(称为 D1)拆分为 D1 内的子补丁 C1,具有相同的以 D1 为中心,使得差异 D1 \ C1 是围绕 C1 的窄框。设上面定义的地图w = w(u), z = v,C1的图像,记为C2。然后在球体 R2 上找到从补丁 D1 到补丁 D2 的变换(映射),其面积与 D1 相同并包含 C2, 你可以将两个子图拼凑成一张图:

w = w(u)
z = v

for [u,v] from C1 ---> [w,z] from C2

w = w_ext(u, v)
z = v

for [u,v] from D1 \ C1 ---> [w,z] from D2 \ C2

问题是如何找到扩展转换w_ext(u)。要使 D2 的面积等于 D1 的面积,您需要选择 w_ext(u) 以便

integra_(D1 \ C1)  sin(w_ext(u)) dw_ext/du(u) du dv = (R1/R2)^2 Area(D1) - Area(C2) ( = the areas on the right are constants )

现在,选择一个合适的函数(如果你愿意,你可以从一个常数开始)f(u),比如一个系数可调的多项式,这样

integra_(D1 \ C1)  f(u) du dv = (R1/R2)^2 Area(D1) - Area(C2) 

e.g.
f(u) = L (constant) such that
integra_(D1 \ C1)  L du dv = (R1/R2)^2 Area(D1) - Area(C2)

i.e.
L = ( (R1/R2)^2 Area(D1) - Area(C2) ) / integra_(D1 \ C1) du dv

然后求解微分方程

sin(w) dw/du = f(u)

e.g.

sin(w) dw/du = L

w(u) = arccos(L*u + a)

但在这种情况下,将此解决方案与前一个解决方案结合起来很重要,因此 w_ext(u) 的初始条件很重要,可能取决于方向 v,即

w_ext(u, v) = arccos(L*u + a(v))

所以有一种比较费力的方法,但是细节很多,也比较复杂。