在python中，在二维零网格上的两点之间插入一行值

Question

TLDR：

在 2d numpy 数组中找到 2 个点后，如何在 0 数组中在它们之间插入一行 1？

上下文：

目前我正在尝试根据二值化医学图像数据（0 和 1）对 3d 数组进行 2d 操作。最终目标是在填充的体素/像素（即第一个和最后一个实例）的起点和终点之间添加一条 1 线。

为此，我使用 SimpleITK 切出单行，然后将其转换为 numpy 数组。在其他示例之后，我编写了返回一组数组的函数，这些数组显示填充 (1) 像素和空 (0) 像素。

从最早和最后一个实例开始，我通过向所有点添加 1 来“填充”这些点，然后将 2 替换为 1。我想要做的是将这条线添加回二维数组，最终，回到 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. 
    """

预期的结果将是一个看起来像原始的数组，但从“填充”像素的开始到结束有一条清晰的线。

Answer 1

因此，感谢 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)