Can the phase of a function be extracted from only its absolute value and its Fourier transform's absolute value?

The short answer is no. Take constant function $f(x) \equiv C \in \mathbb{C}$. Disregarding normalization, we have $\hat{f} = C \delta$ (in the sense of distributions). Clearly, there is no way to recover the phase of $C$ once we take the absolute value on both sides.

To make this a little more explicit, consider a lot easier version of the problem on the group $G = \mathbb{Z} / N\mathbb{Z}$. Its dual is $\hat{G} = G$. If you'll write out the Fourier transform equations (i.e. $\hat{f}(k) = \sum_{n=0}^{N-1} f(n) \exp(-{i k n \over 2 \pi})$), you'll obtain $2N$ real equations for $2N$ coefficients ($N$ Fourier phases and $N$ original phases). The properties of this system of equations are not clear to me, but the case $N=1$ (this is the same as in the first paragraph, but here we don't need to talk about distributions) already shows that the solutions need not be unique.

I hope someone else can provide more information, I'd be also interested to see what conditions on $f$ one needs to assume to get a unique solution. Even for the case $G = \mathbb{Z} / N\mathbb{Z}$ this looks interesting enough.

The right way to ask the question is: given a function $f\in L²(\mathbb{R})$, can $f$ be determined from $|f|$ and $|\widehat{f}|$ up to a multiplicative constant $c$ of modulus $|c|=1$.

This question dates back to Pauli and the answer is no. One can construct counter examples of the form $a\gamma(x-x_0)+b\gamma(x)+c\gamma(x+x_0)$ with $a,b,c$ properly chose ($\gamma(x)=e^{-\pi x²}$ the standard gaussian so that it is ots own Fourier transform). An other construction is as follows:

take $\chi=\mathbf{1}_{[0,1/2]}$ $(a_j)_{j\in\mathbb{Z}}$ a sequence with finite support (to simplify) and $f(x)=\sum_j a_j\chi(x-j)$ so that $\hat f(\xi)=\sum_j a_je^{2i\pi j\xi}\hat\chi(\xi)$.

Now we want to construct a sequence $(b_j)$ such that $|a_j|=|b_j|$ and $\left|\sum_j a_je^{2i\pi j\xi}\right|=\left|\sum_j b_je^{2i\pi j\xi}\right|$. This can be done via a Riesz product: take $\alpha_1,\ldots,\alpha_N$ a finite real sequence, $\varepsilon_1,\ldots,\varepsilon_N$ a finite sequence of $\pm1$ and consider $$ \prod_{k=1}^N (1+i\alpha_j\varepsilon_j\sin 2\pi 3^j\xi)=\sum a_j^{(\varepsilon)}e^{2i\pi j\xi}. $$

Changing a $\varepsilon_j$ from $+1$ to $-1$ conjugates one of the factors on the left hand side, so it does not change the modulus. Now the same happens for the $a_j^{(\varepsilon)}$: each of them is either $0$ or a product of $i\alpha_j\varepsilon_j$ (up to a constant) -- the point is that it is not the sum of products of $i\alpha_j\varepsilon_j$'s, this is why we took the $3^j$ in the sine!.

This is the phase retrieval problem, for which simple iterative, numerical algorithms exist. For an overview, see J. R. Fienup, "Phase retrieval algorithms: a comparison," Appl. Opt. 21, 2758-2769 (1982).