【python】获取三维图形的某个截面

转自：http://blog.csdn.net/huozi07/article/details/50544575

 在数据可视化实践过程中经常需要对三维甚至更高纬度数据进行可视化。由于视线阻挡，人们在看三维物体时并不能观测清楚完全。有时候需要获取三维图形的某个截面来单独分析数据。

[python]  view plain  copy

 

1. # -*- coding: utf-8 -*-  
2. import numpy as np  
3. from mpl_toolkits.mplot3d import Axes3D  
4. import matplotlib as mpl  
5. from scipy import interpolate  
6. import matplotlib.cm as cm  
7. import matplotlib.pyplot as plt  
8.   
9. xmin=-10  
10. ymin=-10  
11. xmax=10  
12. ymax=10  
13. rawStrideNum=50#原始数据x,y的分段数  
14. InterpolationNum=100#使用插值方法获取某个截面数据  获取数据时将数据分割的分段数  
15. #生成原始数据的f value  
16. def fm((x, y)):  
17.     return (np.sin(x) + 0.05 * x ** 2+ np.sin(y) + 0.05 * y ** 2)  
18. #原始数据获取  
19. x = np.linspace(xmin, xmax, rawStrideNum)  
20. y = np.linspace(ymin,ymax,rawStrideNum)  
21. X, Y = np.meshgrid(x, y)#50*50的网格数据  
22. Z = fm((X, Y))  
23.   
24. #开始作图  
25. fig = plt.figure(figsize=(9, 6))  
26. #Draw sub-graph1  
27. ax=plt.subplot(1, 2, 1,projection = '3d')  
28. surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True)  
29. ax.set_xlabel('x')  
30. ax.set_ylabel('y')  
31. ax.set_zlabel('f(x, y)')  
32. plt.colorbar(surf, shrink=0.5, aspect=5)#标注  
33.   
34. #subplot2data preparation  
35. #(1)用插值函数拟合数据  
36. newfun=interpolate.interp2d(X, Y, Z,kind="cubic")#返回的是一个函数  newfun为一个函数  
37. xnew = np.linspace(xmin, xmax, InterpolationNum)#重新划分x,y间隔数为InterpolationNum个  
38. ynew = np.linspace(ymin, ymax, InterpolationNum)  
39. fnew=newfun(xnew,ynew)#得到重新划分x,y间隔后fnew值  fnew为InterpolationNum*InterpolationNum的二维数组  
40.   
41. def getXCrossSectionData(f,xval):  
42.     index=(xval-xmin)*1.0/(xmax-xmin)*InterpolationNum  
43.     print " x index  ",int(index)  
44.     FixedX=np.empty(InterpolationNum)  
45.     FixedX.fill(xval)  
46.     return  FixedX,f[int(index)]  
47. def getYCrossSectionData(f,yval):  
48.     index2=(yval-ymin)*1.0/(ymax-ymin)*InterpolationNum  
49.     print "y index",int(index2)  
50.     FixedY=np.empty(InterpolationNum)  
51.     FixedY.fill(yval)  
52.     return FixedY,f[:,int(index2)]  
53.   
54. FixedX,Z_FixedX=getXCrossSectionData(fnew,-9)#设置想要获取的x截面  
55. FixedY,Z_FixedY=getYCrossSectionData(fnew,5)  
56. #在原图上添加直线，从吻合性判断数据的准确性  
57. ax.plot(FixedX, ynew, Z_FixedX)  
58. ax.plot(xnew, FixedY, Z_FixedY)  
59.   
60. #Draw sub-graph2单独抽取某个截面  
61. ax2=plt.subplot(122, projection = '3d')  
62. ax2.plot(FixedX, ynew, Z_FixedX)  
63. ax2.plot(xnew, FixedY, Z_FixedY)  
64.   
65. plt.show()  

    由上图可知，左图在描绘三维物品时，因为视觉阻挡，并不能很好观察底部数据。通过抽取数据的某些截面，可以单独分析数据。