Filling functions on hypercubes given the function's integrals over each coordinate

Question

$\newcommand\dif{\mathop{}\!\mathrm{d}}$ Suppose $f:[0,1]^2\rightarrow\mathbb{R}_{>0}$ and that: $$\int_0^1{f(x,y)\dif x}=a(y),\hspace{1cm}\int_0^1{f(x,y)\dif y}=b(x),$$ for some known functions $a,b:[0,1]\rightarrow\mathbb{R}_{>0}$.

Obviously, there are many possible $f$ consistent with these two integral equations, but if this is really all we know, one "natural" solution for $f$ is that: $$f(x,y)=m\frac{a(y)}{m}\frac{b(x)}{m},$$ where: $$m:=\int_0^1{\int_0^1{f(x,y)\dif x}\dif y}=\int_0^1{b(x)\dif x}=\int_0^1{a(y)\dif y}.$$

Question 1: Can you formalize the notion of "naturalness" here? E.g., is there some (natural!) constrained minimisation problem with this solution?

Now suppose instead that $f:[0,1]^3\rightarrow\mathbb{R}_{>0}$ and that: $$\int_0^1{f(x,y,z)\dif x}=a(y,z),\hspace{1cm}\int_0^1{f(x,y,z)\dif y}=b(x,z),\hspace{1cm}\int_0^1{f(x,y,z)\dif z}=c(x,y),$$ for some known functions $a,b,c:[0,1]^2\rightarrow\mathbb{R}_{>0}$.

Question 2: What is the equivalent "natural" solution for $f$ now?

Question 3: How does this generalize to dimensions above 3?

---

Edit 1: I had been looking for answers that were a product as in the 2D case. However, it seems that the "natural" answers are in fact linear combinations of products instead.

In the 3D case, let: $$p(x):=\int_0^1{\int_0^1{f(x,y,z)\dif y}\dif z}=\int_0^1{b(x,z)\dif z}=\int_0^1{c(x,y)\dif y},$$ $$q(y):=\int_0^1{\int_0^1{f(x,y,z)\dif x}\dif z}=\int_0^1{a(y,z)\dif z}=\int_0^1{c(x,y)\dif x},$$ $$r(z):=\int_0^1{\int_0^1{f(x,y,z)\dif x}\dif y}=\int_0^1{a(y,z)\dif y}=\int_0^1{b(x,z)\dif x},$$ $$m:=\int_0^1{\int_0^1{\int_0^1{f(x,y,z)\dif x}\dif y}\dif z}=\int_0^1{p(x)\dif x}=\int_0^1{q(y)\dif y}=\int_0^1{r(z)\dif z}.$$ Then one solution is: $$f(x,y,z)=m\frac{p(x)}{m}\frac{a(y,z)}{m}+m\frac{q(y)}{m}\frac{b(x,z)}{m}+m\frac{r(z)}{m}\frac{c(x,y)}{m}-2m\frac{p(x)}{m}\frac{q(y)}{m}\frac{r(z)}{m}.$$

Question 4: Is this positive, despite the minus sign?

---

Edit 2: Some evidence that it may not always be positive. Consider the discrete analogue of this problem, where integrals are replaced by sums. Applying the previous formula to fill in the following array (where 1 indices correspond to sums over indices 2 and 3):

val(:,:,1) =
        8.368       4.0451       4.3229
       3.5973       1.4781       2.1193
       4.7706        2.567       2.2036
val(:,:,2) =
       3.1518      0.82929       2.3225
      0.55961            ?            ?
       2.5922            ?            ?
val(:,:,3) =
       5.2162       3.2158       2.0004
       3.0377            ?            ?
       2.1785            ?            ?

gives:

val(:,:,1) =
        8.368       4.0451       4.3229
       3.5973       1.4781       2.1193
       4.7706        2.567       2.2036
val(:,:,2) =
       3.1518      0.82929       2.3225
      0.55961     -0.12621      0.68582
       2.5922      0.95549       1.6367
val(:,:,3) =
       5.2162       3.2158       2.0004
       3.0377       1.6043       1.4334
       2.1785       1.6115      0.56698

which has a negative entry. Note: This discrete problem can be nested inside the original integral one by considering functions that are constant on $[0,1/2)$ and $[1/2,1]$ with respect to each coordinate.

Question 5: Is there another (more natural?) solution which is always positive?