Как сопоставить значения разных сеток? - PullRequest
Новая
       102

Как сопоставить значения разных сеток?

1 голос
/ 13 июля 2020

Я пытаюсь создать набор изображений цилиндра в оттенках серого, исходя из информации о радиусе, размере и ориентации цилиндра.

Этот цилиндр должен содержаться в кубе c 3D me sh, где каждая точка сетки будет представлять пиксель.

Код будет примерно таким. 

  #Define meshgrid
  x_ = np.linspace(0,Lx,int(Lx/vox))
  y_ = np.linspace(0,Ly,int(Ly/vox))
  z_ = np.linspace(0,Lz,int(Lz/vox))

  X,Y,Z = np.meshgrid(x_,y_,z_,indexing='ij')

Где vox - это размер вокселя, который я ввел. Это определит 3D me sh, Матрица, из которой мы получим стек шкалы серого, будет определена как 

  M = np.zeros((x_.size,y_.size,z_.size),dtype=int)

Теперь нам дан цилиндр, центральная ось которого идет от точки p1 к p2, и имеет радиус r. Итак, на основе следующей ссылки Numpy маска из координат цилиндра , я создаю дополнительную сетку, содержащую все точки, принадлежащие этому цилиндру. 

        #Vector
        v = p2 - p1

        #Normalize vector
        lenght = scipy.linalg.norm(v)
        v = v / lenght

        # make some vector not in the same direction as v
        not_v = np.array([1.0, 0, 0])
        if (v == not_v).all():
              not_v = np.array([0, 1.0, 0])
        # make vector perpendicular to v
        n1 = np.cross(v, not_v)
        # normalize n1
        n1 = n1 / scipy.linalg.norm(n1)
        # make unit vector perpendicular to v and n1
        n2 = np.cross(v, n1)

        #Define gridpoints for the cilinder
        l_ =      np.linspace(0,lenght,100)
        r_ =      np.linspace(0,r,10)
        theeta_ = np.linspace(0,2*np.pi,10)

        #define meshgrid for cilinder
        L,R,Theeta = np.meshgrid(l_,r_,theeta_,indexing='ij')

Тогда я бы получил координаты xyz по этим цилиндрическим координатам. 

#Transform to x, y, z coordinates
Xc, Yc, Zc = [p1[i] + v[i] * L + R * np.sin(Theeta) * n1[i] + r * np.cos(Theeta) * n2[i] for i in [0, 1, 2]]

Итак, теперь вопрос в следующем. Я определил точки сетки, которые составляют данный цилиндр в декартовых координатах. Однако проблема, с которой я сейчас сталкиваюсь, заключается в попытке сопоставить эту сетку для цилиндра с сеткой, определенной для получения стека изображений в градациях серого.

Итак, я хотел бы взять этот X c , Y c, Z c, координаты и посмотреть, каким координатам X, Y, Z они будут соответствовать, и осветить соответствующие пиксели в матрице M. Но я не вижу очевидного способа сделать это.

С уважением,

1 Ответ

1 голос
/ 13 июля 2020

Вы неправильно подходите к задаче, вам следует напрямую найти, какие точки в вашей координатной сетке попадают в ваш цилиндр. Вот как это сделать: 

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

# Problem parameters
Lx, Ly, Lz = 50, 20, 30
vox = 2.0
p1 = np.array([12., 8., 4.])
p2 = np.array([37., 14., 22.])
r = 5.
#Define meshgrid
x_ = np.linspace(0, Lx, int(Lx / vox))
y_ = np.linspace(0, Ly, int(Ly / vox))
z_ = np.linspace(0, Lz, int(Lz / vox))
# Grid coordinates
X, Y, Z = np.meshgrid(x_, y_, z_, indexing='ij')
# Stack into an array
coords = np.stack([X, Y, Z], axis=-1)
# Compute distance from each point to cylinder axis
v = p2 - p1
t = np.dot(coords - p1, v) / np.dot(v, v)
p = p1 + np.expand_dims(t, axis=-1) * v
dist = np.linalg.norm(coords - p, axis=-1)
# Select points within cylinder distance and bounds
mask = (dist <= r) & (0 <= t) & (t <= 1)
print(mask.shape)
# (10, 4, 6)
# Select coordinates in cylinder
cyl_coords = coords[mask]
# Plot
ax = plt.figure().add_subplot(111, projection='3d')
ax.scatter3D(cyl_coords[:, 0], cyl_coords[:, 1], cyl_coords[:, 2], s=3)
# Set plot limits for uniform aspect ratio
ax.set_xlim(0, 50)
ax.set_ylim(-15, 35)
ax.set_zlim(-10, 40)
ax.figure.tight_layout()
ax.figure.show()

Вывод:

Выходное изображение

...