在两个图像之间插值

Question

我正在尝试在 Python 中的两个图像之间进行插值。

图像具有形状 (188, 188)

我希望在这两个图像之间插入图像。假设 Image_1 位于位置 z=0，而 Image_2 位于位置 z=2。我想要位置 z=1 的插值图像。

我相信这个答案（MATLAB）包含类似的问题和解决方案。

使用 MATLAB 在 3D MRI 体积中创建中间切片

我尝试将此代码转换为 Python，如下所示：

from scipy.interpolate import interpn
from scipy.interpolate import griddata

# Construct 3D volume from images
#    arr.shape = (2, 182, 182)
arr = np.r_['0,3', image_1, image_2]

slices,rows,cols = arr.shape

# Construct meshgrids
[X,Y,Z] = np.meshgrid(np.arange(cols), np.arange(rows), np.arange(slices));
[X2,Y2,Z2] = np.meshgrid(np.arange(cols), np.arange(rows), np.arange(slices*2));

# Run n-dim interpolation
Vi = interpn([X,Y,Z], arr, np.array([X1,Y1,Z1]).T)

但是，这会产生错误：

ValueError: The points in dimension 0 must be strictly ascending

我怀疑我没有正确构建我的网格，但是我对这种方法是否正确有点迷茫。

有任何想法吗？

- - - - - 编辑 - - - - - -

找到了一些似乎可以解决此问题的 MATLAB 代码：

在 3d 空间中的两个平面之间进行插值

我试图将其转换为 Python：

from scipy.ndimage.morphology import distance_transform_edt
from scipy.interpolate import interpn

def ndgrid(*args,**kwargs):
    """
    Same as calling ``meshgrid`` with *indexing* = ``'ij'`` (see
    ``meshgrid`` for documentation).
    """
    kwargs['indexing'] = 'ij'
    return np.meshgrid(*args,**kwargs)

def bwperim(bw, n=4):
    """
    perim = bwperim(bw, n=4)
    Find the perimeter of objects in binary images.
    A pixel is part of an object perimeter if its value is one and there
    is at least one zero-valued pixel in its neighborhood.
    By default the neighborhood of a pixel is 4 nearest pixels, but
    if `n` is set to 8 the 8 nearest pixels will be considered.
    Parameters
    ----------
      bw : A black-and-white image
      n : Connectivity. Must be 4 or 8 (default: 8)
    Returns
    -------
      perim : A boolean image

    From Mahotas: http://nullege.com/codes/search/mahotas.bwperim
    """

    if n not in (4,8):
        raise ValueError('mahotas.bwperim: n must be 4 or 8')
    rows,cols = bw.shape

    # Translate image by one pixel in all directions
    north = np.zeros((rows,cols))
    south = np.zeros((rows,cols))
    west = np.zeros((rows,cols))
    east = np.zeros((rows,cols))

    north[:-1,:] = bw[1:,:]
    south[1:,:]  = bw[:-1,:]
    west[:,:-1]  = bw[:,1:]
    east[:,1:]   = bw[:,:-1]
    idx = (north == bw) & \
          (south == bw) & \
          (west  == bw) & \
          (east  == bw)
    if n == 8:
        north_east = np.zeros((rows, cols))
        north_west = np.zeros((rows, cols))
        south_east = np.zeros((rows, cols))
        south_west = np.zeros((rows, cols))
        north_east[:-1, 1:]   = bw[1:, :-1]
        north_west[:-1, :-1] = bw[1:, 1:]
        south_east[1:, 1:]     = bw[:-1, :-1]
        south_west[1:, :-1]   = bw[:-1, 1:]
        idx &= (north_east == bw) & \
               (south_east == bw) & \
               (south_west == bw) & \
               (north_west == bw)
    return ~idx * bw

def signed_bwdist(im):
    '''
    Find perim and return masked image (signed/reversed)
    '''    
    im = -bwdist(bwperim(im))*np.logical_not(im) + bwdist(bwperim(im))*im
    return im


def bwdist(im):
    '''
    Find distance map of image
    '''
    dist_im = distance_transform_edt(1-im)
    return dist_im


def interp_shape(top, bottom, num):
    if num<0 and round(num) == num:
        print("Error: number of slices to be interpolated must be   integer>0")

    top = signed_bwdist(top)
    bottom = signed_bwdist(bottom)

    r, c = top.shape
    t = num+2

    print("Rows - Cols - Slices")
    print(r, c, t)
    print("")

    # rejoin top, bottom into a single array of shape (2, r, c)
    # MATLAB: cat(3,bottom,top)
    top_and_bottom = np.r_['0,3', top, bottom]
    #top_and_bottom = np.rollaxis(top_and_bottom, 0, 3)

    # create ndgrids 
    x,y,z = np.mgrid[0:r, 0:c, 0:t-1]  # existing data
    x1,y1,z1 = np.mgrid[0:r, 0:c, 0:t] # including new slice

    print("Shape x y z:", x.shape, y.shape, z.shape)
    print("Shape x1 y1 z1:", x1.shape, y1.shape, z1.shape)
    print(top_and_bottom.shape, len(x), len(y), len(z)) 

    # Do interpolation
    out = interpn((x,y,z), top_and_bottom, (x1,y1,z1))

    # MATLAB: out = out(:,:,2:end-1)>=0;
    array_lim = out[-1]-1    
    out[out[:,:,2:out] >= 0] = 1 

    return out

我这样称呼它：

new_image = interp_shape(image_1,image_2, 1)

我很确定这是那里的 80% 但我在运行时仍然收到此错误：

ValueError: The points in dimension 0 must be strictly ascending

同样，我可能没有正确构建我的网格。我相信np.mgrid应该产生与 MATLABs 相同的结果ndgrid。

有没有更好的方法来构造ndgrid等价物？

Answer 1

我想通了。或者至少是一种产生理想结果的方法。

基于：在 3d 空间中的两个平面之间进行插值

def signed_bwdist(im):
    '''
    Find perim and return masked image (signed/reversed)
    '''    
    im = -bwdist(bwperim(im))*np.logical_not(im) + bwdist(bwperim(im))*im
    return im

def bwdist(im):
    '''
    Find distance map of image
    '''
    dist_im = distance_transform_edt(1-im)
    return dist_im

def interp_shape(top, bottom, precision):
    '''
    Interpolate between two contours

    Input: top 
            [X,Y] - Image of top contour (mask)
           bottom
            [X,Y] - Image of bottom contour (mask)
           precision
             float  - % between the images to interpolate 
                Ex: num=0.5 - Interpolate the middle image between top and bottom image
    Output: out
            [X,Y] - Interpolated image at num (%) between top and bottom

    '''
    if precision>2:
        print("Error: Precision must be between 0 and 1 (float)")

    top = signed_bwdist(top)
    bottom = signed_bwdist(bottom)

    # row,cols definition
    r, c = top.shape

    # Reverse % indexing
    precision = 1+precision

    # rejoin top, bottom into a single array of shape (2, r, c)
    top_and_bottom = np.stack((top, bottom))

    # create ndgrids 
    points = (np.r_[0, 2], np.arange(r), np.arange(c))
    xi = np.rollaxis(np.mgrid[:r, :c], 0, 3).reshape((r**2, 2))
    xi = np.c_[np.full((r**2),precision), xi]

    # Interpolate for new plane
    out = interpn(points, top_and_bottom, xi)
    out = out.reshape((r, c))

    # Threshold distmap to values above 0
    out = out > 0

    return out


# Run interpolation
out = interp_shape(image_1,image_2, 0.5)

示例输出：

Answer 2

我不知道你的问题的解决方案，但我认为不可能用interpn.

我更正了您尝试的代码，并使用了以下输入图像：

但结果是：

这是更正后的代码：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy import interpolate

n = 8
img1 = np.zeros((n, n))
img2 = np.zeros((n, n))

img1[2:4, 2:4] = 1
img2[4:6, 4:6] = 1

plt.figure()
plt.imshow(img1, cmap=cm.Greys)

plt.figure()
plt.imshow(img2, cmap=cm.Greys)

points = (np.r_[0, 2], np.arange(n), np.arange(n))
values = np.stack((img1, img2))
xi = np.rollaxis(np.mgrid[:n, :n], 0, 3).reshape((n**2, 2))
xi = np.c_[np.ones(n**2), xi]

values_x = interpolate.interpn(points, values, xi, method='linear')
values_x = values_x.reshape((n, n))
print(values_x)

plt.figure()
plt.imshow(values_x, cmap=cm.Greys)
plt.clim((0, 1))

plt.show()

我认为你的代码和我的代码之间的主要区别在于xi. interpn使用起来往往有些混乱，我在一个旧的答案中更详细地解释了它。如果您对我如何指定的机制感到好奇xi，请参阅我的这个答案，解释我所做的事情。

这个结果并不完全令人惊讶，因为interpn只是在两个图像之间进行线性插值：因此一个图像中具有 1 且另一图像中具有 0 的部分简单地变成了 0.5。

在这里，由于一个图像是另一个图像的翻译，很明显，我们想要一个在“中间”翻译的图像。但是如何interpn插入两个一般图像呢？如果你有一个小圆圈和一个大圆圈，那么在它们之间“”应该有一个中间大小的圆圈是否清楚？在狗和猫之间进行插值怎么样？或者狗和建筑物？

我认为你本质上是在尝试“画线”连接两个图像的边缘，然后尝试找出之间的图像。这类似于以半帧采样移动视频。您可能想查看诸如光流之类的东西，它使用向量连接相邻帧。我不知道是否以及有哪些 python 包/实现可用。