python - 在python中,在二维零网格上的两点之间插入一行值
问题描述
TLDR:
在 2d numpy 数组中找到 2 个点后,如何在 0 数组中在它们之间插入一行 1?
语境:
目前我正在尝试从二值化医学图像数据(0 和 1)对 3d 数组进行 2d 操作。最终目标是在填充的体素/像素(即第一个和最后一个实例)的起点和终点之间添加一条 1 线。
为此,我使用 SimpleITK 切出单行,然后将其转换为 numpy 数组。在其他示例之后,我编写了返回一组数组的函数,这些数组显示填充 (1) 像素和空 (0) 像素。
从最早和最后一个实例开始,我通过向所有点添加 1 来“填充”这些点,然后将 2 替换为 1。我想要做的是将这条线添加回 2d 数组,并最终返回到3d 阵列。
我知道我有一个通过 scipy.ndimage.map_coordinates(np.transpose(z), line_array) 返回的坐标列表,但我不知道如何将其应用于原始数组。
在我看来,我可以简单地看到创建一个 0 的二维数组,然后可以简单地将其添加到原始二维数组中。即便如此,我还是不知道如何将线插入 2d 数组,然后插入 3d 数组。首先在 2d 中工作的愿望是因为这些数组可能非常非常大。任何帮助将不胜感激。
import numpy as np
import scipy
import matplotlib.pyplot as plt
import simpleITK as sitk
from timeit import default_timer as timer
#Functions used are shown below
#Read in a CT scan using SimpleITK
ctscan = sitk.ReadImage("example.mhd")
#Extract a slice along the Z axis (Simple ITK uses x, y, z indexing)
z = ctscan[:,:,150:151]
#Convert the slice to a numpy array, which is then z, y, x indexing
z = sitk.GetArrayFromImage(z)
#Drop the 3rd dimension to view in matplotlib
z = z[0, :, :]
plt.imshow(z)
#Get the line and the line array
line, line_array = extract_line(0, 0, z.shape[1], z.shape[0], z)
#Returned from scipy.ndimage.map_coordinates
#1d array of 490 elements
#In [153]: len(line)
#Out[153]: 490
#2d tuple with coordinates
#In [154]: line_array.shape
#Out[154]: (2, 490)
#Get the indices that are filled, the first and last contiguous set
points, start, end = get_line_limits(line)
#Number of arrays that are filled with 1s
#In [158]: len(points)
#Out[158]: 36
#First instance of contiguous 1's
#In [159]: start
#Out[159]: array([35, 36, 37], dtype=int64)
#Last instance of contiguous 1's
#In [160]: end
#Out[160]: array([424, 425], dtype=int64)
#If I am interpreting this then that should be the first point
#In [161]: x1 = line_array[0][35]
#In [162]: y1 = line_array[1][35]
#And this the last
#In [161]: x2 = line_array[0][425]
#In [162]: y2 = line_array[1][425]
#From here I know that I can create a "blank" array, but I don't know how
#To "place" the line...
V = np.zeros((z.shape[0], z.shape[1]))
In [166]: V = np.zeros((z.shape[0], z.shape[1]))
'''
In [167]: V
Out[167]:
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])
'''
#and then
new = z + V
#Hopefully this would then show the line between the two defined points #array.
plt.imshow(new)
"""
Functions. It's likely that the answer is already somewhere in here but I've missed it from sheer idiocy.
"""
def extract_line(x0, y0, x1, y1, z):
"""
Extract a line from a 2d slice.
:param x0: Starting point of line on x axis, usually 0
:param y0: Starting point of line on y axis, usually 0
:param x1: Ending point of line on x axis, usually the len
:param y1: Ending point of line on x axis, usually the len
:param z: The 2d Numpy array
:return: Returns an interpolated line with filled pixels and their coordinates. Taken from https://stackoverflow.com/questions/7878398/how-to-extract-an-arbitrary-line-of-values-from-a-numpy-array
"""
start = timer()
x0, y0 = float(x0), float(y0) # These are in pixel coordinates!!
x1, y1 = float(x1), float(y1)
x_len = abs(x0 - x1)
y_len = abs(y0 - y1)
line_length = int(np.sqrt((x_len**2) + (y_len**2))) #Length of line
x, y = np.linspace(x0, x1, line_length), np.linspace(y0, y1, line_length)
line_array = np.vstack((x, y))
# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(np.transpose(z), line_array)
'''
#Uncomment this is you want to view the results.
#Also I understand matplotlib clearly does this somehow, any way to
#Simple return the array?
fig, axes = plt.subplots(nrows=2)
axes[0].imshow(z)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')
axes[1].plot(zi)
plt.show()
'''
stop = timer()
print(abs(start - stop))
return zi, line_array
def consecutive(data, stepsize=1):
"""
Function to find consective elements in an array.
:param data: numpy array.
:param stepsize: how many values between elements before splitting the array.
:return: Returns an array broken along the step size.
"""
consecutive = np.split(data, np.where(np.diff(data) != stepsize)[0]+1)
return consecutive
def get_line_limits(line):
"""
Function to find the first and last instance of the filed pixels as determined by extract_line.
:param line: numpy array returned from extract line.
:return: Returns the indeces of the filled (i.e. 1) pixels, along with the first and last set.
"""
start = timer()
line_points = np.where(line==1)
#Call the consecutive function to get the contiguous points.
line_points = consecutive(line_points[0])
line_start = line_points[0]
line_end = line_points[(len(line_points) - 1)]
stop = timer()
print("It took {} seconds to get limits.".format(abs(start-stop)))
return line_points, line_start, line_end
def place_line(x0, y0, x1, y1, z):
"""
Mystical function that doesn't exist yet because I can't seem to work it out.
"""
预期的结果将是一个看起来像原始的数组,但从“填充”像素的开始到结束有一条清晰的线。
解决方案
因此,感谢 Ardweaden 建议我找到一种简单的方法来陈述我的目标,以及 Mad Physicist 为我指出了 Bresenham 线算法的方向。在那之后,我找到了一些弥补差距的 python 实现。两个版本的工作代码如下(如果有人有更快的方法,我会将其标记为已接受的答案):
#Starting from the line limits obtained in the original questions code
x0 = np.ceil(line_array[0][start[0]]) #first filled element on x
y0 = np.ceil(line_array[1][start[0]]) #first filled element on y
x1 = np.floor(line_array[0][end[(len(end) - 1)]]) #last filled element on x
y1 = np.floor(line_array[1][end[(len(end) - 1)]]) #last filled element on y
#The skimage implementation of Bresenham's line algorithm
# https://scikit-image.org/docs/dev/api/skimage.draw.html#skimage.draw.line
from skimage.draw import line
#Make a 2d array of 0s that is the same size as the original 2d array z
V = np.zeros((z.shape[0], z.shape[1]), dtype=np.uint8)
#This flips the y and x from how I think of it...
rr, cc = line(int(y0), int(x0), int(y1), int(x1))
#Cast the line coordinates as 1s
V[rr, cc] = 1
#Then add the new array to the old array (although could just write the line on z).
new_z = V + z
#Use numpy where to change the 2 values back to 1 values.
new_z[np.where(new_z==2)] = 1
#Show the final image
plt.imshow(new_z)
# This implementation is a short script from https://github.com/encukou/bresenham and can just be
# brought into the existing code base without importing an external library.
# The code can be modified, but as is doesn't have a return so place it into a list
new_line = list(bresenham(int(y0), int(x0), int(y1), int(x1)))
#Use list comprehension to get the arrays
x_array = np.array([item[0] for item in new_line])
y_array = np.array([item[1] for item in new_line])
#Then cast them as 1's
V[x_array, y_array] = 1
#Combine the original and the new array
new_z = V + z
#Recast the 2's to 1's
new_z[np.where(new_z==2)] = 1
#And view
plt.imshow(new_z)
推荐阅读
- python - 带有 get_object 的 UserPassesTestMixin 的 test_func 不适用于 ListView
- c# - 如何从使用实体框架的外键链接的多个表中获取所有数据?
- node.js - 如何将 DynamoDB Streams 对象映射到 Javascript 对象?
- r - 更改格子图例的字体
- java - Spring MVC:在同一网络服务器中的 WAR 之间共享 bean
- php - 登录问题。我正在使用 php 版本 7.2
- python - 如何将 tweepy api 调用转换为异步
- javascript - 在 ListView 中显示 Array 的每一项
- c# - 无法使用标准序列化时如何将 C# 对象传输到 C++
- react-redux - Redux-saga - yield all 不是在 root saga 中调用我所有的 saga