使用imshow或其他函数创建交叉图像

Question

我有两个三维的地面穿透雷达数据数组。每个数组其实就是一系列随时间变化的二维图像，时间是在第三个维度上增加的。我想创建一个三维图形，这个图形可以同时显示每个数组中的一幅二维图像。

我基本上是在尝试制作一个围栏图。这种类型的图可以在以下网站找到一些例子： http://www.geogiga.com/images/products/seismapper_3d_seismic_color.gif http://www.usna.edu/Users/oceano/pguth/website/so461web/seismic_refl/fence.png

我通常使用imshow来单独显示这些二维图像以进行分析。但是，我研究了一下imshow的功能，发现它似乎不能与三维坐标轴一起使用。有没有什么方法可以解决这个问题？或者有没有其他的绘图函数可以实现imshow的功能，并且可以与三维坐标轴结合使用？

Answer 1

用matplotlib是无法做到这一点的。@DrV的回答只是个近似值。matplotlib实际上并不会显示原始图像的每一个像素，而是显示了一张经过重新缩放的图像。rstride和cstride可以帮助你指定图像缩放的方式，但输出的结果并不是完全和原图一样的。

Answer 2

可能还有更好的方法，但至少你可以先创建一个平面网格并给它上色：

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

# create a 21 x 21 vertex mesh
xx, yy = np.meshgrid(np.linspace(0,1,21), np.linspace(0,1,21))

# create some dummy data (20 x 20) for the image
data = np.random.random((20, 20))

# create vertices for a rotated mesh (3D rotation matrix)
X = np.sqrt(1./3) * xx + np.sqrt(1./3) * yy
Y = -np.sqrt(1./3) * xx + np.sqrt(1./3) * yy
Z = np.sqrt(1./3) * xx - np.sqrt(1./3) * yy

# create the figure
fig = plt.figure()

# show the reference image
ax1 = fig.add_subplot(121)
ax1.imshow(data, cmap=plt.cm.BrBG, interpolation='nearest', origin='lower', extent=[0,1,0,1])

# show the 3D rotated projection
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=plt.cm.BrBG(data), shade=False)

这样就会生成：

（请注意，我在旋转矩阵上没有太小心，你需要自己创建投影。使用一个真正的旋转矩阵可能是个好主意。）

另外要注意的是，围栏的柱子和围栏之间有个小问题，也就是说，网格的顶点比补丁的数量多了一个。

---

上面的方法在处理高分辨率图像时效率不高，甚至可能没什么用。另一种选择是使用支持仿射图像变换的后端。不幸的是，这样你就得自己计算变换。虽然这并不是特别困难，但还是有点麻烦，而且你无法得到一个可以旋转的真实3D图像等等。

关于这种方法，可以参考 http://matplotlib.org/examples/api/demo_affine_image.html

另外，你可以使用OpenCV和它的 cv2.warpAffine 函数，在用 imshow 显示图像之前对图像进行变形。如果你把周围填充成透明颜色，就可以叠加图像，得到一个看起来像你示例图的效果。

---

为了让你了解 plot_surface 的可能性，我尝试把Lena的图像变形为一个半圆柱：

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

# create a 513 x 513 vertex mesh
xx, yy = np.meshgrid(np.linspace(0,1,513), np.linspace(0,1,513))

# create vertices for a rotated mesh (3D rotation matrix)
theta = np.pi*xx
X = np.cos(theta)
Y = np.sin(theta)
Z = yy

# create the figure
fig = plt.figure()

# show the 3D rotated projection
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=plt.imread('/tmp/lena.jpg')/255., shade=False)

她确实变形得很好，但对图像的所有操作都比较慢：

Answer 3

如果你愿意尝试其他的绘图库（也就是说，不用matplotlib），那么可以考虑mayavi或tvtk这两个库，虽然学习起来可能会有点难。根据我看到的，最接近你想要的效果的是标量切割平面，具体可以参考这个链接：http://wiki.scipy.org/Cookbook/MayaVi/Examples

大部分的文档可以在这里找到：http://docs.enthought.com/mayavi/mayavi/index.html