The standard variational derivation of stellar matter structure in the Wigner-Seitz approximation is generalized to the finite temperature situation where a wide distribution of different nuclear species can coexist in the same density and proton fraction condition, possibly out of $\beta$-equilibrium. The same theoretical formalism is shown to describe on one side the single-nucleus approximation (SNA), currently used in most core collapse supernova simulations, and on the other side the nuclear statistical equilibrium (NSE) approach, routinely employed in r- and p-process explosive nucleosynthesis problems. In particular we show that in-medium effects have to be accounted for in NSE to have a theoretical consistency between the zero and finite temperature modeling. The bulk part of these in-medium effects is analytically calculated and shown to be different from a van der Waals excluded volume term. This unified formalism allows controlling quantitatively the deviations from the SNA in the different thermodynamic conditions, as well as having a NSE model which is reliable at any arbitrarily low value of the temperature, with potential applications for neutron star cooling and accretion problems. We present different illustrative results with several mass models and effective interactions, showing the importance of accounting for the nuclear species distribution even at temperatures lower than 1 MeV.