I have two equations of images
$$g_1=f_1+h_2*f_2$$
and
$$g_2=f_2+h_1*f_1$$
where $*$ is convolution and $g_1,g_2,h_1,h_2$ are known.
How do I find $f_1$ and $f_2$?
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.