This booklet addresses your questions in Chapter 0.9 and specifically in equation (6). The Moyal equation is just the von Neumann evolution equation for the density matrix in a different, phase space, representation. The von Neumann equation amounts to the Schroedinger equation, i.e. they specify time development of functions, with absolutely no information on what such functions are. If the initial condition inputs are physical wavefunctions $\psi(x,t_0)$ or Wigner distributions (WF) $f(x,p,t_0)$ with all the requisite requirements for physicality (normalization, reality for the WF, eqn (6), etc), the evolution output is a physical function at time t. If the inputs are are garbage, the outputs are most likely garbage.
You are thus probably asking about the requisite conditions for a phase space function to be physical at a specific moment, e.g. the origin of time. For a real normalized phase-space function to be of the form
$$
f(x,p)=\frac{1}{2\pi}\int\! dy~\psi^* \left (x-\frac{\hbar}{2}
y \right )~e^{-iyp} ~ \psi \left (x+\frac{\hbar}{2} y \right ),
$$
so bona-fide physical and normalized, you need to ensure its cross-spectral density, its Fourier transform, ``left-right" factorizes,
$$
\tilde{f}(x,y)=\int dp ~e^{ipy} f(x,p) ~ = ~
g^{*}_L (x-\hbar y/2) ~g_R (x+\hbar y/2)~.
$$
That is,
$$
{\partial^2 ~~~\ln \tilde{f} \qquad \qquad \phantom{a} \over
\partial(x-\hbar y/2)~\partial(x+\hbar y/2)}
=0 ~, \tag{6}
$$
so that, for real $f$, ~ $g_L=g_R$.
This all holds for pure states; for mixed states, you have to parse and reorganize the above argument.