some question about B-spline kernel matrix

[MangoWAY]

Hi! I'm studying the source code. I find the method

def bspline_kernel_3d(sigma=[1, 1, 1], order=2, asTensor=False, dtype=th.float32, device='cpu'):
    kernel_ones = th.ones(1, 1, *sigma)
    kernel = kernel_ones
    padding = np.array(sigma) - 1

    for i in range(1, order + 1):
        kernel = F.conv3d(kernel, kernel_ones, padding=(padding).tolist())/(sigma[0]*sigma[1]*sigma[2])
	
    if asTensor:
        return kernel[0, 0, ...].to(dtype=dtype, device=device)
    else:
        return kernel[0, 0, ...].numpy()

I know we can use the control points to control the displacements. I've also read about other bspline implementations. The airlab bspline method looks like very simple. So I want to know how does it work, the meaning of sigma and why the for loop can get the bspline kernel. Thanks! :)

[cherise215]

Not sure if my answer is 100 per cent correct or not. But in WIKI it says: 'Fast b-spline interpolation on a uniform sample domain can be done by iterative mean-filtering'. https://en.wikipedia.org/wiki/B-spline. In this case, the initial kernel can be view as a mean-filter, which is then interactively convolved to achieve higher-order interpolation.

[qiuhuaqi]

Just adding some details to @cherise215's answer to explain this implementation:

• n-th order B-spline function can be computed by recursive convolution with 0-th order B-spline function n times: 0-th order B-spline function is a rectangle window, which is kernel_one in the code. And sigma is effectively the "width" of this rectangle.

• Since this rectangle function is symmetrical, convolution between kernel and kernel_one is the same as cross-correlation, which can be implemented using F.convNd() function. By setting padding=sigma-1, only the over-lapping part of the convolution is computed, which is also referred to as the "support" (non-zero part) of B-spline function.

• You should also expect the size of the kernel getting bigger as the order increases. The kernel size should increase by sigma-1 in each convolution, so output kernel size should be sigma + order*(sigma-1). You can work out the sigma you want from your desired kernel size using this formula.

As an example, 2D quadratic (order=2) B-spline kernel with sigma=[6, 6] (kernel size 16x16) looks like this

• This kernel is then used in https://github.com/airlab-unibas/airlab/blob/80c9d487c012892c395d63c6d937a67303c321d1/airlab/transformation/pairwise.py#L624 to interpolate control points parameters to acquire dense displacements.

This implementation allows us to use an arbitrary order of B-spline function, which is pretty neat.

[qiuhuaqi]

A small problem in my opinon is that the transformation model in this implementation doesn't seem to allow flexible combinations of kernel size and control point spacing (stride=sigma in transposed conv), which are both locked by sigma.