How to efficiently implement a sub-pixel convolution layer in a CNN for image super-resolution

I am playing around using convolution neural networks for image super-resolution, i.e. up-scaling of images. In literature I often come across the concept of sub-pixel convolution as one of the layers in CNNs for super-resolution. If you want to upscale by a factor r, the idea is to merge r^2 feature maps by interleaving them into a single image. Below is an example of 4 feature maps of 3x3 being merged into a single 6x6 image (i.e. scaling factor of 2). This can be extended to bigger scaling factors.

I have been looking at available neural network layers in MMA 12.1, but I have not figured out how to implement a sub-pixel convolution with the available layers. Maybe the use of a deconvolution layer with a clever choice of a fixed kernel and strides might do the trick but I can't see how yet.

Did anybody in the MMA community already try to implement a sub-pixel convolution efficiently, or maybe point me in some direction how to implement it using the available NN layers.

You are right that the sub-pixel convolution can be achieved with deconvolution layers with fixed weights

net = NetGraph[{
"part1" -> PartLayer[1],
"part2" -> PartLayer[2],
"part3" -> PartLayer[3],
"part4" -> PartLayer[4],
"dcov1" ->
DeconvolutionLayer["Weights" -> {{{{1, 0}, {0, 0}}}},
"Biases" -> None, "Stride" -> 2],
"dcov2" ->
DeconvolutionLayer["Weights" -> {{{{0, 1}, {0, 0}}}},
"Biases" -> None, "Stride" -> 2],
"dcov3" ->
DeconvolutionLayer["Weights" -> {{{{0, 0}, {1, 0}}}},
"Biases" -> None, "Stride" -> 2],
"dcov4" ->
DeconvolutionLayer["Weights" -> {{{{0, 0}, {0, 1}}}},
"Biases" -> None, "Stride" -> 2],
"sum" -> TotalLayer[]
}, {"part1" -> "dcov1" -> "sum", "part2" -> "dcov2" -> "sum",
"part3" -> "dcov3" -> "sum", "part4" -> "dcov4" -> "sum"}]


data = {{{{2, 1, 5}, {4, 0, 9}, {3, 0, 0}}}, {{{4, 3, 0}, {0, 1,
7}, {2, 8, 4}}}, {{{5, 8, 9}, {1, 3, 2}, {0, 1, 6}}}, {{{0, 4,
1}, {7, 2, 3}, {4, 5, 3}}}};

net[data][[1]] // MatrixForm