This book has many programs for conjugate-pair transformations. More transformations are described in pseudocode and more multidimensional transformations are found in my imaging book. Operators can be cascaded with one another yielding a greater multitude of conjugate-pair operators. From this multitude the missing data program shown used the convolution subroutine contran() but another operator from the multitude could have been chosen. The convolution operator was a good choice because convolution amounts to the family of ordinary differential equations--a big family of applications. What other operators could lead to useful results? This is food for thought.
Obviously the same missing data program could be used for two-dimensional data if we replaced the one-dimensional convolution subroutine by two-dimensional convolution. Simple filters could be interpreted as partial differential equations such as for plate bending or for a tent over poles. Large filters would be more costly, but they might be represented as a separable approximation. Theoretically the CG method could require many iterations--as many as there are missing values. This creates the potential for burdensome computational costs but convergence is often much faster than the theoretical limit. If you have a convergence speed problem you might get some guidance from the literature on solving elliptic equations.
What the missing data program does is provide missing values in data space by spreading out residuals in another space. Let us think about other transformations and spaces.