How to solve a two variable equation of images?
If we add hats to denote the Fourier transforms of those functions, convolution turns into multiplication, and addition stays as addition:

$$\hat g_1=\hat f_1+\hat h_2\times\hat f_2$$
$$\hat g_2=\hat f_2+\hat h_1\times\hat f_1$$

If we feed that to Wolfram alpha, sans hats:

solve({g_1=f_1+h_2*f_2, g_2=f_2+h_1*f_1}, {f_1, f_2})


the output is, with hats added back:

$$\hat f_1 = \frac{\hat g_2 \hat h_2-\hat g_1}{\hat h_1 \hat h_2-1} \text{ and } \hat f_2 = \frac{\hat g_1 \hat h_1-\hat g_2}{\hat h_1 \hat h_2-1} \text{ and } \hat h_1 \hat h_2-1\ne0$$

Frequency domain division equals time domain deconvolution. The frequency domain $1$ represents a time domain identity convolution kernel of a single 1 at $(0, 0)$ and all other values zero. Reading from the above, the solution is valid only if $h_2$ convolved by $h_2$ minus the identity kernel gives an invertible filter, in other words one without frequency domain zeros.