14

実験データの 3D 表現を構築して、膜の変形を追跡したいと考えています。実験的には、コーナー ノードのみがわかっています。ただし、構造全体の変形をプロットしたいので、膜を補間してその素晴らしいカラーマップを有効にしたいのです。周りを検索することで、次のコードでほぼそれに近づきました。

import numpy
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.interpolate import griddata

x=numpy.array([0, 0, 1, 1])
y=numpy.array([0.5, 0.75, 1, 0.5])
z=numpy.array([0, 0.5, 1,0])

fig = plt.figure()
ax = Axes3D(fig)
verts = [zip(x, y, z)]
PC = Poly3DCollection(verts)
ax.add_collection3d(PC)

xi = numpy.linspace(x.min(),x.max(),20)
yi = numpy.linspace(y.min(),y.max(),20)
zi = griddata((x,y),z, (xi[None,:], yi[:,None]), method='linear')
xig, yig = numpy.meshgrid(xi, -yi)
ax.plot_surface(xig, yig, zi, rstride=1, cstride=1,  linewidth=0,cmap=plt.cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))
plt.show()

次のプロットを取得します。

ここに画像の説明を入力

青いポリゴンは、そのコーナー ノードによって認識されるサーフェスであり、カラーマップしたいものです。カラーマップされたサーフェスは、これまでのところ私の最高の結果です。しかし、私を悩ませている表面の上部近くに黒いポリゴンがあります。表面がメッシュグリッドに適合していないため、4番目のコーナーがナンになっている可能性があると思います。

これらの黒い三角形を回避するための回避策はありますか、またはそのコーナー ノードだけが認識しているサーフェスをカラーマッピングするより良い方法はありますか?

編集:これは、次のコマンドを使用して最初のコメントで指定された三角測量ソリューションの図です。

triang = tri.Triangulation(x, y)
ax.plot_trisurf(x, y, z, triangles=triang.triangles, cmap=cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))

ここに画像の説明を入力

4

2 に答える 2

6

確かにplot_trisurf、このタスクには完璧なはずです! さらに、 を使用して、より小さな三角形tri.UniformTriRefinerを取得できます。Triangulation

import numpy
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import tri, cm

x = numpy.array([0, 0, 1, 1])
y = numpy.array([0.5, 0.75, 1, 0.5])
z = numpy.array([0, 0.5, 1, 0])

triang = tri.Triangulation(x, y)
refiner = tri.UniformTriRefiner(triang)
new, new_z = refiner.refine_field(z, subdiv=4)

norm = plt.Normalize(vmax=abs(y).max(), vmin=-abs(y).max())
kwargs = dict(triangles=new.triangles, cmap=cm.jet, norm=norm, linewidth=0.2)

fig = plt.figure()
ax = Axes3D(fig)
pt = ax.plot_trisurf(new.x, new.y, new_z, **kwargs)
plt.show()

次の画像になります。

ここに画像の説明を入力

Triangular Grid Refinement は最近追加されたばかりmatplotlibなので、使用するにはバージョン 1.3 が必要です。ただし、バージョン 1.2 で行き詰まる場合は、行とすべてのメソッドをコメントアウトすれば、Github からソースを直接使用することもできるはずです。次に、メソッドと使用を使用して、新しい対応する Z 値を補間する必要があります。import matplotlib.tri.triinterpolaterefine_fieldrefine_triangulationgriddata

編集:上記のコードは、3 次補間を使用して新しい三角形の Z 値を決定しますが、線形補間の場合は、次の行を代用または追加できます。

interpolator = tri.LinearTriInterpolator(triang, z)
new, new_z = refiner.refine_field(z, interpolator, subdiv=4)

または、で補間を行うにはscipy.interpolate.griddata：

from scipy.interpolate import griddata

new = refiner.refine_triangulation(subdiv = 4)
new_z = griddata((x,y),z, (new.x, new.y), method='linear')
于 2013-11-15T23:51:23.677 に答える
4

問題は、matplotlib でサーフェスの補間シェーディングを行う方法、つまり Matlab のshading('interp')機能に相当する方法に要約されます。簡単に言えば、できません。ネイティブでサポートされていないため、これまでに提示されたソリューションが目指しているのは、手動で行うことです。

数年前、Matlab にも不満を感じていたときに、この道をたどりshading('interp')ました。各四辺形の 4 つのコーナーの色を補間するだけで機能します。つまり、隣接する四辺形では色のグラデーションの方向が異なる可能性があります。私が望んでいたのは、各カラー バンドが正確に z 軸上の 2 つの明確に定義された値の間にあり、隣接するセル間に視覚的な切れ目がないことでした。

三角測量に取り組むことは間違いなく正しい考えです。しかし、単純にグリッドを改良し、隣接する三角形の色が視覚的に区別できなくなるポイントに到達することを期待する代わりに (アーティファクトが最初に表示されるポイントに到達することなく)、私のアプローチは、三角形分割の輪郭バンドを計算し、それらを 3D でプロットすることでした。 .

これを最初に実装したとき、matplotlib は三角形分割の輪郭をサポートしていませんでした。今では経由で行い_tri.TriContourGeneratorます。これが、抽出されたポリゴン頂点の z 値も提供していれば、完了です。残念ながら、これらは Python レベルではアクセスできないため、 と の出力を比較して再構築する必要がありますcreate_filled_contours()。create_contours()これは次のコードで行われます。

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from matplotlib import _tri, tri, cm

def contour_bands_3d(x, y, z, nbands=20):
    # obtain the contouring engine on a triangulation
    TRI = tri.Triangulation(x, y)
    C = _tri.TriContourGenerator(TRI.get_cpp_triangulation(), z)

    # define the band breaks
    brks = np.linspace(z.min(), z.max(), nbands+1)

    # the contour lines
    lines = [C.create_contour(b) for b in brks]

    # the contour bands
    bands = [C.create_filled_contour(brks[i], brks[i+1]) for i in xrange(nbands)]

    # compare the x, y vertices of each band with the x, y vertices of the upper
    # contour line; if matching, z = z1, otherwise z = z0 (see text for caveats)
    eps = 1e-6
    verts = []
    for i in xrange(nbands):
        b = bands[i][0]
        l = lines[i+1][0]
        z0, z1 = brks[i:i+2]
        zi = np.array([z1 if (np.abs(bb - l) < eps).all(1).any() else z0 for bb in b])
        verts.append(np.c_[b, zi[:,None]])
    return brks, verts

x = np.array([0, 0, 1, 1])
y = np.array([0.5, 0.75, 1, 0.5])
z = np.array([0, 0.5, 1,0])

fig = plt.figure()
ax = Axes3D(fig)
verts = [zip(x, y, z)]
PC = Poly3DCollection(verts)
ax.add_collection3d(PC)

# calculate the 3d contour bands
brks, verts = contour_bands_3d(x, -y, z)

cmap = cm.get_cmap('jet')
norm = plt.Normalize(vmax=abs(y).max(), vmin=-abs(y).max())

PC = Poly3DCollection(verts, cmap=cmap, norm=norm, edgecolors='none')
PC.set_array(brks[:-1])
ax.add_collection(PC)
ax.set_ylim((-1, 1))
plt.show()

結果は次のとおりです。

縞模様のある膜

z 値の再構成は完全には正しくないことに注意してください。これは、ax、y 頂点が実際に元のデータセットの一部であるかどうかも確認する必要があるためです。その場合、元の z 値を取得する必要があります。ただし、コンター アルゴリズムの C++ コードを変更して z 値を追跡する方がはるかに簡単です。これは小さな変更ですが、Python ですべてのケースをカバーしようとするのは悪夢にほかなりません。

効率に関しては、まあ、Python レベルでグラフィックス カードの仕事をしようとしているので、ひどいものになるでしょう。しかし、それはすべて同じですmplot3d。パフォーマンスの実装が必要な場合は、 VTKBandedContourFilter()をお勧めします。これは非常に高速に動作し、Python からも使用できます。

于 2013-11-17T02:51:58.803 に答える