Let $N$ be a natural number. Consider the following set of matrices whose entries are non-negative integers: $$X_N:=\left\{(c_{ij})_{i,j=1}^4\in M_4(\mathbb{Z}_{\geq 0})\bigg| \sum_j c_{1j} = \sum_i c_{2j}, \sum_i c_{i1}=\sum_i c_{i2}, \sum_{i,j}c_{ij} = N\right\}.$$ Write $f(N) = |X_N|$.

Question: Is there an algorithm / theorem which gives an explicit formula for $f(N)$ as a function of $N$? It should be polynomial. Without the constraints about the sum of the elements of the first two columns is equal and the sum of the elements of the first two rows is equal, this will just be $\binom{N+15}{15}$ which is a polynomial of degree 15. I am looking for a similar formula for the number of elements of $X_N$.