import numpy as np
def view_field(height_map, viewer_height, viewer_pos, resolution=360):
height_map = np.asarray(height_map)
h, w = height_map.shape
vi, vj = viewer_pos
# Find locations higher than viewer
m = height_map > viewer_height
# Find angle of each pixel relative to viewer
ii, jj = np.ogrid[-vi:h - vi, -vj:w - vj]
a = np.arctan2(ii, jj)
# Distance of each pixel to viewer
d2 = np.square(ii) + np.square(jj)
d = np.sqrt(d2)
# Find angle range "behind" each pixel
pix_size = 0.5
ad = np.arccos(d / np.sqrt(d2 + np.square(pix_size)))
# Minimum and maximum angle encompassed by each pixel
amin = a - ad
amax = a + ad
# Define angle "bins"
ar = np.linspace(-np.pi, np.pi, resolution + 1)
# Find the bin corresponding to each pixel
b = np.digitize(a, ar) % resolution
bmin = np.digitize(amin, ar) % resolution
bmax = np.digitize(amax, ar) % resolution
# Find the closest distance to a high pixel for each angle bin
angdist = np.full_like(ar, np.inf)
np.minimum.at(angdist, bmin[m], d[m])
np.minimum.at(angdist, bmax[m], d[m])
# Visibility is true if the pixel distance is less than the
# visibility distance for its angle bin
return d <= angdist[b]
下面是它的运行方式:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# A wavy function for testing
def ackley(x, y):
return (-20 * np.exp(-0.2 * np.sqrt(0.5 * (np.square(x) + np.square(y))))
- np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y)))
+ np.e + 20)
# Make input data
x, y = np.ogrid[-4:4:100j, -4:4:100j]
height_map = ackley(x, y)
viewer_height = 7.5
# Viewer coordinates are in "pixel space"
viewer_pos = (50, 70)
# Compute visibility
f = view_field(height_map, viewer_height, viewer_pos)
# Plot height map and result
plt.figure(figsize=(8, 4))
ax = plt.subplot(121, projection='3d')
ax.plot_surface(np.arange(100)[:, np.newaxis], np.arange(100),
height_map, color='b', alpha=0.6)
ax.scatter3D(*viewer_pos, viewer_height, s=5, c='r', label='Viewer')
plt.title('Height map')
plt.legend()
plt.subplot(122)
plt.imshow(f)
plt.title('Visibility map')
plt.tight_layout()
我提出了这个函数,它不是很完美,但是做的和你说的类似。它有一个{{CD1}}参数,它表示你想考虑的圆周中有多少个“角度步”。我不确定这是解决这个问题的“正确”方法,而且我觉得应该有一个更聪明的方法来做一些操作,但无论如何
下面是它的运行方式:
结果是:
相关问题 更多 >
编程相关推荐