Suppose we have several interacting particles in pure state $\left|\psi\right>$. For each of particles we can extract density matrices via
$$\rho_i(x_i,x_i^\prime)=\int \left<X_i,x_i\middle|\psi\right>\left<\psi\middle|X_i,x_i^\prime\right>dX_i,$$
where $x_i$ is coordinates of particle $i$ and $X_i$ is coordinates of the rest of system.
Can the original wave function $\left<X\middle|\psi\right>$ of the whole system be restored from all the $\rho_i$?
No, it is not possible in general, because the particles are almost surely to be entangled (since they are interacting, as mentioned). The reason is that each $\rho_i$ is nothing but the reduced density operator of the $i$ th particle. Indeed, in the form you have written, $$\rho_i(x_i,x_i') = \langle x_i | \rho_i | x_i' \rangle = \langle x_i | \mathsf{Tr}_{X_i \neq x_i}(|\psi \rangle \langle \psi |) | x_i' \rangle ,$$ where the $X_i$ go over the rest of the system. Even if all the $\rho_i$ are known, one cannot find $\rho \equiv |\psi \rangle \langle \psi |$, the state of the system, because there can exist correlations among the variables of different particles that cannot permit one to write the following: $$ \rho = \rho_1 \otimes \rho_2 \otimes .. \rho_i \otimes .. \ \ \ \ \mathsf{(Incorrect)} $$ This is just the continuous variable version of entanglement. As an example with qubits, consider the Bell state $$ |\Psi \rangle = \frac{|01\rangle + |10\rangle}{\sqrt{2}} $$ and trace out one of the qubits. This yields $$ \rho_1 = \rho_2 = \frac{|0\rangle\langle 0|+|1\rangle\langle 1|}{2}, $$ But even a state with only classical correlations, $$ \rho = \frac{|00\rangle \langle00| +|11\rangle\langle11|}{2} $$ would give the same reduced density operators, this state being much different from the entangled one considered. This shows that there is more information in the full state $\rho$ than in the partial knowledge of each of the subsystems $\rho_i$.
As Abhinav states, in general, it is not possible. Assuming low interaction you can however start with a Weyl transform
$$\begin{align*} P_i(x,p) &= \frac{1}{\pi\hbar}\int_{-\infty}^\infty \langle x+y| \hat{\rho}_i |x-y \rangle e^{-2ipy/\hbar}\,dy &\Bigg|\quad \hat\rho_i =\iint du\,dv\,\rho_i(u,v)|u\rangle\langle v| \\ &= \frac1{\pi\hbar}\int_{-\infty}^\infty \rho_i(x+y,x-y) e^{-2ipy/\hbar}\,dy \end{align*}$$
Then use $|\psi_i(x)|^2 = \int_{-\infty}^\infty P_i(x,p)\,dp$ and $|\phi_i(p)|^2=\int_{-\infty}^\infty P_i(x,y)\,dx$ to obtain the intensity of the individual wave functions in coordinate- and momentum-space. Finally, reconstruct the phase e.g. via Gerchberg-Saxton. Obviously this will fail horribly once sufficient interaction/mixture occured...
