Given
$\mathrm R \in \mathbb R^{m \times n}$, whose $m$ rows are the differences of neighboring rows of $\mathrm X \in \mathbb R^{m \times n}$
$\mathrm C \in \mathbb R^{m \times n}$, whose $n$ columns are the differences of neighboring columns of $\mathrm X \in \mathbb R^{m \times n}$
we would like to determine the unknown matrix $\mathrm X \in \mathbb R^{m \times n}$. Do note that we consider the 1st and $m$-th rows neighbors, and that we also consider the 1st and $n$-th columns neighbors.
The subtraction matrix
We define the $d \times d$ circulant matrix
$$\mathrm S_d := \begin{bmatrix} -1 & 1 & 0 & \dots & 0 & 0\\ 0 & -1 & 1 & \dots & 0 & 0\\ 0 & 0 & -1 & \dots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \dots & -1 & 1\\ 1 & 0 & 0 & \dots & 0 & -1\end{bmatrix}$$
Let
$$w_d := \exp \left( i \frac{2 \pi}{d} \right)$$
From results on circulant matrices [0], we know that the spectrum of $\mathrm S_d$ is
$$\{ w_d^0 - 1, w_d^1 - 1, \dots, w_d^{d-1} - 1 \} = \{ 0, w_d - 1, \dots, w_d^{d-1} - 1 \}$$
In other words, the $d$ eigenvalues of $\mathrm S_d$ are on the unit circle centered at $-1 + i 0$. Hence, we can conclude that $\mbox{rank} (\mathrm S_d) = d-1$. Since $\mathrm S_d 1_d = 0_d$, the (right) null space of $\mathrm S_d$ is $\{ \gamma 1_d \mid \gamma \in \mathbb R \}$.
A linear system
We have two linear matrix equations in $\mathrm X$
$$\mathrm S_m \mathrm X = \mathrm R \qquad\qquad\qquad \mathrm X \mathrm S_n^{\top} = \mathrm C$$
which is a linear system of $2 m n$ scalar equations in $m n$ unknowns. Note that
- if $\bar{\mathrm X}$ is a particular solution of $\mathrm S_m \mathrm X = \mathrm R$, then so is $\bar{\mathrm X} + \gamma \, 1_m \mathrm v^{\top}$, where $\gamma \in \mathbb R$ and $\mathrm v \in \mathbb R^n$, as
$$\mathrm S_m \left( \bar{\mathrm X} + \gamma \, 1_m \mathrm v^{\top} \right) = \underbrace{\mathrm S_m \bar{\mathrm X}}_{= \mathrm R} + \gamma \, \underbrace{\mathrm S_m 1_m}_{= 0_m} \mathrm v^{\top} = \mathrm R + \mathrm O_{m \times n} = \mathrm R$$
- if $\bar{\mathrm X}$ is a particular solution of $\mathrm X \mathrm S_n^{\top} = \mathrm C$, then so is $\bar{\mathrm X} + \gamma \, \mathrm u 1_n^{\top}$, where $\gamma \in \mathbb R$ and $\mathrm u \in \mathbb R^m$, as
$$\left( \bar{\mathrm X} + \gamma \, \mathrm u 1_n^{\top} \right) \mathrm S_n^{\top} = \underbrace{\bar{\mathrm X} \mathrm S_n^{\top}}_{= \mathrm C} + \gamma \, \mathrm u \underbrace{1_n^{\top} \mathrm S_n^{\top}}_{= \mathrm 0_n^{\top}} = \mathrm C + \mathrm O_{m \times n} = \mathrm C$$
Thus, if $\bar{\mathrm X}$ is a particular solution of both $\mathrm S_m \mathrm X = \mathrm R$ and $\mathrm X \mathrm S_n^{\top} = \mathrm C$, then so is
$$\bar{\mathrm X} + \gamma \, 1_m 1_n^{\top}$$
where $\gamma \in \mathbb R$. Since there are infinitely many solutions, more information is needed to reconstruct $\mathrm X$.
Suppose now that we are given not only matrices $\mathrm R$ and $\mathrm C$, but also $s \in \mathbb R$, which denotes the sum of the $m n$ entries of $\mathrm X$. Appending the equation $1_m^{\top} \mathrm X 1_n = s$ to the linear system
$$\mathrm S_m \mathrm X = \mathrm R \qquad\qquad\qquad \mathrm X \mathrm S_n^{\top} = \mathrm C \qquad\qquad\qquad 1_m^{\top} \mathrm X 1_n = s$$
and vectorizing, we obtain a linear system of $2 m n + 1$ equations in $m n$ unknowns
$$\begin{bmatrix} \mathrm I_n \otimes \mathrm S_m\\ \mathrm S_n \otimes \mathrm I_m\\ 1_{mn}^{\top}\end{bmatrix} \mbox{vec} (\mathrm X) = \begin{bmatrix} \mbox{vec} (\mathrm R)\\ \mbox{vec} (\mathrm C)\\ s\end{bmatrix}$$
Minimizing the Frobenius norm
Let us look for the solution that has the least Frobenius norm. We have the quadratic program (QP)
$$\begin{array}{ll} \text{minimize} & \| \mathrm X \|_{\text{F}}^2\\ \text{subject to} & \mathrm S_m \mathrm X = \mathrm R\\ & \mathrm X \mathrm S_n^{\top} = \mathrm C\end{array}$$
We define the Lagrangian
$$\mathcal L (\mathrm X, \Lambda_1, \Lambda_2) := \frac 12 \| \mathrm X \|_{\text{F}}^2 - \langle \Lambda_1, \mathrm S_m \mathrm X - \mathrm R \rangle - \langle \Lambda_2, \mathrm X \mathrm S_n^{\top} - \mathrm C \rangle$$
Taking the partial derivative with respect to $\mathrm X$ and finding where it vanishes, we obtain
$$\mathrm X = \mathrm S_m^{\top} \Lambda_1 + \Lambda_2 \mathrm S_n$$
Left-multiplying by $\mathrm S_m$ and right-multiplying by $\mathrm S_n^{\top}$, we obtain two coupled linear matrix equations in the Lagrange multipliers $\Lambda_1$ and $\Lambda_2$
$$\begin{array}{rl} \mathrm S_m \mathrm S_m^{\top} \Lambda_1 + \mathrm S_m \Lambda_2 \mathrm S_n &= \mathrm R\\ \mathrm S_m^{\top} \Lambda_1 \mathrm S_n^{\top} + \Lambda_2 \mathrm S_n \mathrm S_n^{\top} &= \mathrm C\end{array}$$
Vectorizing, we obtain a linear system of $2 m n$ equations in $2 m n$ unknowns
$$\begin{bmatrix} \mathrm I_n \otimes \mathrm S_m \mathrm S_m^{\top} & \mathrm S_n^{\top} \otimes \mathrm S_m\\ \mathrm S_n \otimes \mathrm S_m^{\top} & \mathrm S_n \mathrm S_n^{\top} \otimes \mathrm I_m\end{bmatrix} \begin{bmatrix} \mbox{vec} (\Lambda_1)\\ \mbox{vec} (\Lambda_2)\end{bmatrix} = \begin{bmatrix} \mbox{vec} (\mathrm R)\\ \mbox{vec} (\mathrm C)\end{bmatrix}$$
Solving this linear system and unvectorizing the solutions, we obtain the optimal Lagrange multipliers $\Lambda_1^*$ and $\Lambda_2^*$. The least-norm solution is then given by
$$\mathrm X_{\text{LN}} := \mathrm S_m^{\top} \Lambda_1^* + \Lambda_2^* \mathrm S_n$$
Reference
[0] Robert M. Gray, Toeplitz and Circulant Matrices: A Review.
