使用坐标变化和 scipy 插值旋转 2D 图像
Rotating a 2D image using change of coordinates and scipy interpolation
我正在尝试旋转我的图像，但好像我的画框根本没有旋转。
以下是我的代码步骤：
1 - 创建倾斜圆盘的图像。
2 - 在 x 和 y 上应用坐标变化以获得新的坐标系。
3 - 在这些坐标处进行二维插值。
问题：使用余弦函数缩小x轴的倾斜部分可以正常工作，但框架的旋转部分根本不起作用。看来我的骨架缩得更厉害了
import numpy as np
import matplotlib.pyplot as plt
from numpy import cos as cos
from numpy import sin as sin
import scipy
from scipy.interpolate import interp1d
def coord_change(qu, qv, inclination_tmp, ang_tmp):
ang = ang_tmp*np.pi/180
inclination = inclination_tmp*np.pi/180
new_qu = (qu*cos(ang)+qv*sin(ang))*cos(inclination)
new_qv = -qu*sin(ang)+qv*cos(ang)
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def imaging_model(rho,I_profile, R_image, **kwargs):
method = kwargs.get('method', 'linear')
xv = np.linspace(-R_image, R_image, len(rho)*2, endpoint=False) # Taille de la fenetre
interpol_index = interp1d(rho, I_profile,kind=method)
X, Y = np.meshgrid(xv, xv)
profilegrid2 = np.zeros(X.shape, float)
current_radius = np.sqrt(X**2 + Y**2)
cond=np.logical_and(current_radius<=max(rho),current_radius>=min(rho)) # Min et max des données
profilegrid2[cond] = interpol_index(current_radius[cond])
return xv, profilegrid2
champ = 80 # mas
diam = 10 # mas
nb_points = 100
rho = np.linspace(0,champ/2,nb_points)
inclin = 20
angle = 0
I_profile = np.zeros(nb_points)
I_profile[rho<=diam/2] = 1
xv, image = imaging_model(rho, I_profile, champ/2)
##############################
x_mod, y_mod, x_tot_mod = coord_change(xv, xv, inclin, angle)
f_image = scipy.interpolate.interp2d(xv, xv, image)
image_mod = f_image(x_mod,y_mod)
plt.figure()
plt.imshow(image_mod,extent=[min(xv),max(xv),min(xv),max(xv)])
小提示：我不想用python里的旋转功能，想自己动手
更新：
我编写了相同的代码，但将倾斜和旋转过程分为两部分以便更好地理解并更好地检查问题可能来自何处：
import numpy as np
import matplotlib.pyplot as plt
from numpy import cos as cos
from numpy import sin as sin
import scipy
from scipy.interpolate import interp1d
def coord_rotation(qu, qv, ang_tmp):
ang = ang_tmp*np.pi/180
new_qu = qu*cos(ang)+qv*sin(ang)
new_qv = -qu*sin(ang)+qv*cos(ang)
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def coord_inclination(qu, qv, inclination_tmp):
inclination = inclination_tmp*np.pi/180
new_qu = qu*cos(inclination)
new_qv = qv
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def imaging_model(rho,I_profile, R_image, **kwargs):
method = kwargs.get('method', 'linear')
xv = np.linspace(-R_image, R_image, len(rho)*2, endpoint=False) # Taille de la fenetre
interpol_index = interp1d(rho, I_profile,kind=method)
X, Y = np.meshgrid(xv, xv)
profilegrid2 = np.zeros(X.shape, float)
current_radius = np.sqrt(X**2 + Y**2)
cond=np.logical_and(current_radius<=max(rho),current_radius>=min(rho)) # Min et max des données
profilegrid2[cond] = interpol_index(current_radius[cond])
return xv, profilegrid2
champ = 80 # mas
diam = 10 # mas
nb_points = 100
rho = np.linspace(0,champ/2,nb_points)
inclin = 30
angle = 40
I_profile = np.zeros(nb_points)
I_profile[rho<=diam/2] = 1
xv, image = imaging_model(rho, I_profile, champ/2)
##############################
x_inc, y_inc, x_tot_inc = coord_inclination(xv, xv, inclin)
f_inc = scipy.interpolate.interp2d(xv, xv, image)
image_inc = f_inc(x_inc,y_inc)
x_rot, y_rot, x_tot_rot = coord_rotation(xv, xv, angle)
f_rot = scipy.interpolate.interp2d(xv, xv, image_inc)
image_rot = f_rot(x_rot,x_rot)
plt.figure()
plt.imshow(image, extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('Original')
plt.figure()
plt.imshow(image_inc, extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('After inclination')
plt.figure()
plt.imshow(image_rot , extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('After inclination + rotation')
正常图片：OK
倾斜30°后：OK
倾斜 30° 和旋转 40° 后：问题
这里有几个问题。在示例中，我稍微降低了分辨率和转换角度，以便更容易发现问题：
champ = 80 # mas
diam = 40 # mas
nb_points = 20
rho = np.linspace(0,champ/2,nb_points)
inclin = 50
angle = 30
问题 1：网格与 coordinate-aligned 网格
xv 在你的代码中是一个一维数组，数据位置是 xv 的每个元素与任何其他元素的任意组合——一个规则的矩形网格，与坐标轴。
倾斜变换时，x_inc和y_inc也是一维数组，数据移动到x_inc的任意元素与y_inc的任意元素的每一个组合。这可以正常工作，因为转换只会影响图像的 y 轴。这意味着对于 y 的每个值，x 坐标仍然是常数，并且对于 x 的每个值，y 坐标是常数。
然而，旋转打破了这一规则。这需要您使用 full-factorial 数据网格。
这意味着您无需将 xv, xv 传递给变换函数，而是需要为其提供一个完整的阶乘网格，例如创建通过 np.meshgrid。下面显示了您的代码实际在做什么，以及它如何解释数据（将分辨率降低到 20 以使数据点可辨别）：
fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')
ax1.scatter(x_inc, y_inc, marker='.', c='b', s=2)
x_inc_full, y_inc_full = np.meshgrid(x_inc, y_inc)
ax1.scatter(x_inc_full, y_inc_full, marker='.', c='b', s=0.1)
ax1.scatter(x_rot, y_rot, marker='.', c='r', s=2)
x_rot_full, y_rot_full = np.meshgrid(x_rot, y_rot)
ax1.scatter(x_rot_full, y_rot_full, marker='.', c='r', s=0.1)
结果：
您可以看到代码采用对角线输入线和 stretches/rotates 它，然后 re-interprets 它作为网格。
如果我们改为在全网格上执行变换，事情开始看起来更合理：
# again, but doing the transformations on the mesh
fig2, ax2 = plt.subplots()
ax2.set_aspect('equal')
x_orig, y_orig = np.meshgrid(xv, xv)
x_incf, y_incf, xtot_inc = coord_inclination(x_orig, y_orig, inclin)
ax2.scatter(x_incf, y_incf, marker='.', c='b', s=0.2)
x_rotf, y_rotf, xtot_rot = coord_rotation(x_orig, y_orig, angle)
ax2.scatter(x_rotf, y_rotf, marker='.', c='r', s=0.2)
结果：
现在，好多了！
除了...
问题 2：转换并不像您认为的那样
如果您查看上图，您会注意到 x-axis 已被压缩，而不是 y 轴，但在您的倾斜图中圆比高更宽。
那是因为您将图形插值到新网格上，然后在规则网格上绘制插值图形。
这意味着您所看到的实际上是您可能打算应用的变换的逆向变换。倾斜变换压缩了 x 轴，但生成的图像沿 x 轴拉伸了相同的量。这在您的原始图片中可见，但随着我使用的更新参数变得明显：
解决方案是反向应用变换：将图像解释为在变换后的网格上（这就是您希望它看起来的样子），然后将其插值到常规网格上以进行显示：
interp_incf = scipy.interpolate.interp2d(x_incf, y_incf, image)
image_incf = interp_incf(xv, xv)
除了...
问题 3：scipy 不喜欢在非 axis-parallel
的网格上或从网格进行插值
我发现进行涉及任意网格的插值的唯一方法是 运行（呃...）循环我们想要插值到的所有点坐标：
# inclined
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_incf, y_incf)]\
).reshape(x_incf.shape)
但是，上面的问题仍然存在问题编号 2（转换的方式是错误的）。但是，我们无法从 插入 一个 arbitrarily-deformed 网格。这意味着转换本身需要反转：
def inv_coord_rotation(qu, qv, ang_tmp):
ang = ang_tmp*np.pi/180
inv_qu = qu*cos(-ang)+qv*sin(-ang)
inv_qv = -qu*sin(-ang)+qv*cos(-ang)
new_q = np.sqrt((inv_qu)**2 + (inv_qv)**2)
return inv_qu, inv_qv, new_q
def inv_coord_inclination(qu, qv, inclination_tmp):
inclination = inclination_tmp*np.pi/180
inv_qu = qu/cos(inclination)
inv_qv = qv
inv_q = np.sqrt((inv_qu)**2 + (inv_qv)**2)
return inv_qu, inv_qv, inv_q
...现在我们终于可以做我们来这里要做的事情了：
# and now for the pictures, with interpolation applied in correct direction:
fig3, axes = plt.subplots(ncols=3)
# interpolator for the original image:
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
# inclined
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_incfi, y_incfi)]\
).reshape(x_incfi.shape)
axes[0].imshow(image_incf, extent=[min(xv),max(xv),min(xv),max(xv)])
# should be identical to the original inclined version, image_inc
# rotated version of the inclined image
interp_incf = scipy.interpolate.interp2d(xv, xv, image_incf)
x_rotfi, y_rotfi, xtot_rotfi = inv_coord_rotation(x_orig, y_orig, angle)
image_rotf = np.array([[interp_incf(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_rotfi, y_rotfi)]\
).reshape(x_rotfi.shape)
axes[1].imshow(image_rotf, extent=[min(xv),max(xv),min(xv),max(xv)])
然而
...我只是尝试将这两个变换结合起来，以避免插入 already-interpolated 图像，这并不简单，因为它需要您在旋转后应用“倾斜”变换 space，然后事情变得很奇怪。
可能可以使用实现一个函数，该函数采用 non-inverted 转换后的像素坐标并从那里计算出逆变换，但实际上，此时我们只是在处理 scipy.interpolate ，类似于我自己的 odyssey with it. I think the general recommendation if you're dealing with images is to use Pillow，但这超出了我自己的经验。
我正在尝试旋转我的图像，但好像我的画框根本没有旋转。
以下是我的代码步骤：
1 - 创建倾斜圆盘的图像。
2 - 在 x 和 y 上应用坐标变化以获得新的坐标系。
3 - 在这些坐标处进行二维插值。
问题：使用余弦函数缩小x轴的倾斜部分可以正常工作，但框架的旋转部分根本不起作用。看来我的骨架缩得更厉害了
import numpy as np
import matplotlib.pyplot as plt
from numpy import cos as cos
from numpy import sin as sin
import scipy
from scipy.interpolate import interp1d
def coord_change(qu, qv, inclination_tmp, ang_tmp):
ang = ang_tmp*np.pi/180
inclination = inclination_tmp*np.pi/180
new_qu = (qu*cos(ang)+qv*sin(ang))*cos(inclination)
new_qv = -qu*sin(ang)+qv*cos(ang)
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def imaging_model(rho,I_profile, R_image, **kwargs):
method = kwargs.get('method', 'linear')
xv = np.linspace(-R_image, R_image, len(rho)*2, endpoint=False) # Taille de la fenetre
interpol_index = interp1d(rho, I_profile,kind=method)
X, Y = np.meshgrid(xv, xv)
profilegrid2 = np.zeros(X.shape, float)
current_radius = np.sqrt(X**2 + Y**2)
cond=np.logical_and(current_radius<=max(rho),current_radius>=min(rho)) # Min et max des données
profilegrid2[cond] = interpol_index(current_radius[cond])
return xv, profilegrid2
champ = 80 # mas
diam = 10 # mas
nb_points = 100
rho = np.linspace(0,champ/2,nb_points)
inclin = 20
angle = 0
I_profile = np.zeros(nb_points)
I_profile[rho<=diam/2] = 1
xv, image = imaging_model(rho, I_profile, champ/2)
##############################
x_mod, y_mod, x_tot_mod = coord_change(xv, xv, inclin, angle)
f_image = scipy.interpolate.interp2d(xv, xv, image)
image_mod = f_image(x_mod,y_mod)
plt.figure()
plt.imshow(image_mod,extent=[min(xv),max(xv),min(xv),max(xv)])
小提示：我不想用python里的旋转功能，想自己动手
更新：
我编写了相同的代码，但将倾斜和旋转过程分为两部分以便更好地理解并更好地检查问题可能来自何处：
import numpy as np
import matplotlib.pyplot as plt
from numpy import cos as cos
from numpy import sin as sin
import scipy
from scipy.interpolate import interp1d
def coord_rotation(qu, qv, ang_tmp):
ang = ang_tmp*np.pi/180
new_qu = qu*cos(ang)+qv*sin(ang)
new_qv = -qu*sin(ang)+qv*cos(ang)
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def coord_inclination(qu, qv, inclination_tmp):
inclination = inclination_tmp*np.pi/180
new_qu = qu*cos(inclination)
new_qv = qv
new_q = np.sqrt((new_qu)**2 + (new_qv)**2)
return new_qu, new_qv, new_q
def imaging_model(rho,I_profile, R_image, **kwargs):
method = kwargs.get('method', 'linear')
xv = np.linspace(-R_image, R_image, len(rho)*2, endpoint=False) # Taille de la fenetre
interpol_index = interp1d(rho, I_profile,kind=method)
X, Y = np.meshgrid(xv, xv)
profilegrid2 = np.zeros(X.shape, float)
current_radius = np.sqrt(X**2 + Y**2)
cond=np.logical_and(current_radius<=max(rho),current_radius>=min(rho)) # Min et max des données
profilegrid2[cond] = interpol_index(current_radius[cond])
return xv, profilegrid2
champ = 80 # mas
diam = 10 # mas
nb_points = 100
rho = np.linspace(0,champ/2,nb_points)
inclin = 30
angle = 40
I_profile = np.zeros(nb_points)
I_profile[rho<=diam/2] = 1
xv, image = imaging_model(rho, I_profile, champ/2)
##############################
x_inc, y_inc, x_tot_inc = coord_inclination(xv, xv, inclin)
f_inc = scipy.interpolate.interp2d(xv, xv, image)
image_inc = f_inc(x_inc,y_inc)
x_rot, y_rot, x_tot_rot = coord_rotation(xv, xv, angle)
f_rot = scipy.interpolate.interp2d(xv, xv, image_inc)
image_rot = f_rot(x_rot,x_rot)
plt.figure()
plt.imshow(image, extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('Original')
plt.figure()
plt.imshow(image_inc, extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('After inclination')
plt.figure()
plt.imshow(image_rot , extent=[min(xv),max(xv),min(xv),max(xv)])
plt.title('After inclination + rotation')
正常图片：OK
倾斜30°后：OK
倾斜 30° 和旋转 40° 后：问题
这里有几个问题。在示例中，我稍微降低了分辨率和转换角度，以便更容易发现问题：
champ = 80 # mas
diam = 40 # mas
nb_points = 20
rho = np.linspace(0,champ/2,nb_points)
inclin = 50
angle = 30
问题 1：网格与 coordinate-aligned 网格
xv 在你的代码中是一个一维数组，数据位置是 xv 的每个元素与任何其他元素的任意组合——一个规则的矩形网格，与坐标轴。
倾斜变换时，x_inc和y_inc也是一维数组，数据移动到x_inc的任意元素与y_inc的任意元素的每一个组合。这可以正常工作，因为转换只会影响图像的 y 轴。这意味着对于 y 的每个值，x 坐标仍然是常数，并且对于 x 的每个值，y 坐标是常数。
然而，旋转打破了这一规则。这需要您使用 full-factorial 数据网格。
这意味着您无需将 xv, xv 传递给变换函数，而是需要为其提供一个完整的阶乘网格，例如创建通过 np.meshgrid。下面显示了您的代码实际在做什么，以及它如何解释数据（将分辨率降低到 20 以使数据点可辨别）：
fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')
ax1.scatter(x_inc, y_inc, marker='.', c='b', s=2)
x_inc_full, y_inc_full = np.meshgrid(x_inc, y_inc)
ax1.scatter(x_inc_full, y_inc_full, marker='.', c='b', s=0.1)
ax1.scatter(x_rot, y_rot, marker='.', c='r', s=2)
x_rot_full, y_rot_full = np.meshgrid(x_rot, y_rot)
ax1.scatter(x_rot_full, y_rot_full, marker='.', c='r', s=0.1)
结果：
您可以看到代码采用对角线输入线和 stretches/rotates 它，然后 re-interprets 它作为网格。
如果我们改为在全网格上执行变换，事情开始看起来更合理：
# again, but doing the transformations on the mesh
fig2, ax2 = plt.subplots()
ax2.set_aspect('equal')
x_orig, y_orig = np.meshgrid(xv, xv)
x_incf, y_incf, xtot_inc = coord_inclination(x_orig, y_orig, inclin)
ax2.scatter(x_incf, y_incf, marker='.', c='b', s=0.2)
x_rotf, y_rotf, xtot_rot = coord_rotation(x_orig, y_orig, angle)
ax2.scatter(x_rotf, y_rotf, marker='.', c='r', s=0.2)
结果：
现在，好多了！
除了...
问题 2：转换并不像您认为的那样
如果您查看上图，您会注意到 x-axis 已被压缩，而不是 y 轴，但在您的倾斜图中圆比高更宽。
那是因为您将图形插值到新网格上，然后在规则网格上绘制插值图形。
这意味着您所看到的实际上是您可能打算应用的变换的逆向变换。倾斜变换压缩了 x 轴，但生成的图像沿 x 轴拉伸了相同的量。这在您的原始图片中可见，但随着我使用的更新参数变得明显：
解决方案是反向应用变换：将图像解释为在变换后的网格上（这就是您希望它看起来的样子），然后将其插值到常规网格上以进行显示：
interp_incf = scipy.interpolate.interp2d(x_incf, y_incf, image)
image_incf = interp_incf(xv, xv)
除了...
问题 3：scipy 不喜欢在非 axis-parallel
的网格上或从网格进行插值我发现进行涉及任意网格的插值的唯一方法是 运行（呃...）循环我们想要插值到的所有点坐标：
# inclined
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_incf, y_incf)]\
).reshape(x_incf.shape)
但是，上面的问题仍然存在问题编号 2（转换的方式是错误的）。但是，我们无法从 插入 一个 arbitrarily-deformed 网格。这意味着转换本身需要反转：
def inv_coord_rotation(qu, qv, ang_tmp):
ang = ang_tmp*np.pi/180
inv_qu = qu*cos(-ang)+qv*sin(-ang)
inv_qv = -qu*sin(-ang)+qv*cos(-ang)
new_q = np.sqrt((inv_qu)**2 + (inv_qv)**2)
return inv_qu, inv_qv, new_q
def inv_coord_inclination(qu, qv, inclination_tmp):
inclination = inclination_tmp*np.pi/180
inv_qu = qu/cos(inclination)
inv_qv = qv
inv_q = np.sqrt((inv_qu)**2 + (inv_qv)**2)
return inv_qu, inv_qv, inv_q
...现在我们终于可以做我们来这里要做的事情了：
# and now for the pictures, with interpolation applied in correct direction:
fig3, axes = plt.subplots(ncols=3)
# interpolator for the original image:
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
# inclined
interp_orig = scipy.interpolate.interp2d(xv, xv, image)
image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_incfi, y_incfi)]\
).reshape(x_incfi.shape)
axes[0].imshow(image_incf, extent=[min(xv),max(xv),min(xv),max(xv)])
# should be identical to the original inclined version, image_inc
# rotated version of the inclined image
interp_incf = scipy.interpolate.interp2d(xv, xv, image_incf)
x_rotfi, y_rotfi, xtot_rotfi = inv_coord_rotation(x_orig, y_orig, angle)
image_rotf = np.array([[interp_incf(x, y) for x, y in zip(xline, yline)]\
for xline, yline in zip(x_rotfi, y_rotfi)]\
).reshape(x_rotfi.shape)
axes[1].imshow(image_rotf, extent=[min(xv),max(xv),min(xv),max(xv)])
然而
...我只是尝试将这两个变换结合起来，以避免插入 already-interpolated 图像，这并不简单，因为它需要您在旋转后应用“倾斜”变换 space，然后事情变得很奇怪。
可能可以使用实现一个函数，该函数采用 non-inverted 转换后的像素坐标并从那里计算出逆变换，但实际上，此时我们只是在处理 scipy.interpolate ，类似于我自己的 odyssey with it. I think the general recommendation if you're dealing with images is to use Pillow，但这超出了我自己的经验。