当面の問題に対処するために、現在のコードには 2 つの問題があります。
- とは異なり
Sphere、ParametricEllipsoidVTK が構築する方法が原因で多くの迷走ポイントがあるため、最終結果は技術的には閉じたサーフェスではありません。これは、pyvista で修正する予定の既知の問題です。clean(tolerance=1e-5)今のところ、楕円体を呼び出すことでこれを修正できます。これにより、これらのvtkMathエラーがポップアップ表示されなくなります。
- 道路のくぼみから通りを差し引くと、を呼び出すたびに平面の法線が逆に
boolean_differenceなります。これにより、平面の交互の側面に隆起ができますが、これはおそらく望んでいるものではありません。street_mesh.flip_normals()これは、ループの最後で呼び出すことで修正できます。
これら 2 つの変更により、非常にゆっくりではありますが、コードが実行されます (理由はよくわかりませんが、内部でブール演算がどのように機能するかもわかりません)。
from random import choice
import numpy as np
import pyvista as pv
def add_pothole():
street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=100, j_resolution=100)
pothole_template = pv.ParametricEllipsoid(10, 5, 3).clean(tolerance=1e-5).triangulate(
points = street_mesh.points
for i in range(4):
pothole_mesh = pothole_template.copy()
pothole_mesh.translate(choice(points))
# "add" the pothole to the street
street_mesh = pothole_mesh.boolean_difference(street_mesh.triangulate())
street_mesh.flip_normals()
street_mesh.plot()
add_pothole()
z(x, y)ただし、根本的な問題に対処するために、道路が最初からサーフェスである場合、最初からこれを行う必要はありません。パーリン ノイズを考慮すると、おそらくパーリン ノイズからサンプリングされたスカラーを使用して平面から開始し、呼び出しwarp_by_scalar()てノイズ変調を実装します。楕円形のバンプを同じスカラーに追加し、最後にワープを 1 ステップで行うことができます (複数のパーリン ノイズ サンプルを重ね合わせるのとまったく同じ方法です)。
このためにz(x, y)は、楕円体のパラメータ化を計算する必要があります。これは簡単なことではありませんが、それほど難しいことでもありません。R原点を中心とした半径を持つ球の方程式が
x^2 + y^2 + z^2= R^2.
で割った場合、これはより優れています（そして標準的です）R^2：
(x/R)^2 + (y/R)^2 + (z/R)^2 = 1.
楕円体を取得する方法は、各軸を線形に変換し、それぞれの半径を変更することです (したがって、すべてではなくなりますR)。楕円体の暗黙の方程式は次のとおりです。
(x/a)^2 + (y/b)^2 + (z/c)^2 = 1.
楕円体の自明でない中心が必要な場合は、各座標を変換する必要があります。
(x - x0)^2/a^2 + (y - y0)^2/b^2 + (z - z0)^2/c^2 = 1.
これは、3 次元の楕円体の一般的な形式であり、その軸がデカルト座標系の軸に向いていると仮定しています。
これは素晴らしいことです。
(z - z0)^2/c^2 = 1 - (x - x0)^2/a^2 - (y - y0)^2/b^2
(z - z0)^2 = c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2
z = z0 +- sqrt(c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2)
プラスとマイナスの部分は、楕円体がz(x, y)関数ではないという事実に対応しています。これは、各(x, y)ポイントに対して2 つの可能な値があるためです。しかし、上面と底面の 2 つの関数から楕円体を作成できます。
最後に問題に戻ると、ランダムな(x0, y0, z0)点を選択して、上記の楕円体の底面 ( になるもの) を選択することができますz = z0 - sqrt(...)。心配しなければならない唯一のことは、(x, y)楕円体がない点では、虚数の平方根が得られるということです。したがって、最初に楕円体の内側にあるポイントを除外する必要があります。これを行う最も簡単な方法は、とにかく平方根を計算し、NaN を破棄することです。
import numpy as np
import pyvista as pv
# denser plane for sharper ellipsoid
street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=1000, j_resolution=1000)
street_mesh['Normals'] *= -1 # make normals point up for warping later
# add Perlin noise for example
freq = [0.03, 0.03, 0.03]
noise = pv.perlin_noise(5, freq, [0, 0, 0])
street_mesh['height'] = [noise.EvaluateFunction(point) for point in street_mesh.points]
# the relevant part: ellipsoid function
x0, y0, z0 = street_mesh.points[street_mesh.points.size//6, :] # or whatever random point you want
x, y, _ = street_mesh.points.T # two 1d arrays of coordinates
a, b, c = 10, 5, 3 # semi-axes
ellipsoid_fun = z0 - np.sqrt(c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2) # RuntimeWarning
keep_inds = ~np.isnan(ellipsoid_fun)
street_mesh['height'][keep_inds] += ellipsoid_fun[keep_inds]
street_mesh.set_active_scalars('height')
# warp by 'height' and do a quick plot
street_mesh.warp_by_scalar().plot()
上記は、虚数の平方根に関する警告を発します (データ内の NaN につながります)。これが気になる場合は、明示的に沈黙させることができます。または、LBYL に従って、平方根を計算する前に自分で値を確認することもできます。
arg = c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2
keep_inds = arg >= 0
ellipsoid = z0 - np.sqrt(arg[keep_inds])
street_mesh['height'][keep_inds] = ellipsoid
よりシャープな楕円体のために増加した平面密度を使用した結果は次のとおりです。
(最初は、numpy ブロードキャストを使用してすべてのくぼみを一度に計算できると提案しましたが、これを破る NaN フィルタリングのために、実際には簡単ではなく、おそらく可能でさえありません。)