Here's what is going on (this is legacy code, and it is tested):

When an array is passed with `multichannel=True`, we ask the user to place the channels in the final dimension. Recall this entire algorithm requires and expects an array of shape `(M, N, Z, P)` with `M, N, Z` being spatial dimensions and `P` being channels. This early parsing is intended to coerce the input into that format.

When `multichannel=True` and `data.ndim == 3`, what we actually have is a 2-D multichannel array. That third dimension is not spatial, and must not be treated as such. We do need to record the current shape, though, so the end result can be returned in the expected way.

So this code block records `dims`, which will be used for reshaping at the end, and then adds a necessary singleton dimension. After which that singleton is placed in the correct third position in the dimension order, as a `Z` axis placeholder.

Lack of C-contiguous arrays do not greatly affect this algorithm, as it solves the linear system without requiring Cython.