Let $ABCD$ be some tetrahedron with incenter $O$, and let $E$ be the incenter of the tetrahedron $OBCD$.
Draw a small sphere around $O$ and select $B'C'D'$ such that this sphere is inscribed in $AB'C'D'$. This is always possible, by selecting three tangent planes through $A$ that enclose the sphere, and intersecting them with a tangent plane behind the sphere. Let $E'$ be the incenter of $OB'C'D'$. It is clear that things can be chosen such that $E\ne E'$, e.g. by making the small sphere small enough that $E$ is outside $OB'C'D'$.
Now, by assumption we have
$$ f(O)=f(A)f(B)f(C)f(D) \qquad f(O)f(B)f(C)f(D)=f(E) $$
and by multiplying these equations and canceling common factors (since by assumption $f(X)\ne 0$) we get
$$ f(O)^2 = f(A)f(E) $$
Similarly,
$$ f(O)^2 = f(A)f(E') $$
Therefore, canceling again, $f(E)=f(E')$.
However, by appropriate similarity transforms, we can put the entire configuration such that $E$ and $E'$ are anywhere we please, and this calculation still shows that $f(E)=f(E')$.
So $f$ must be a constant function.
The equation $f(O)=f(A)f(B)f(C)f(D)$ implies that the common value of $f(X)$ must be a cube root of $1$, and since the values of $f$ are specified to be real, $f(X)=1$ is the only solution.
If the values are allowed to be complex, we get two additional solutions with $f(X)=e^{\pm 2\pi/3 i}$.
The problem specified that $f(X)\ne 0$ everywhere, but relaxing this gives us only one additional solution, the constant function $f(X)=0$. If there's any point $A$ with $f(A)=0$, then for an arbitrary $O\ne A$ we can find $B'C'D'$ as above such that $O$ is the incenter of $AB'C'D'$. Then $f(O)=f(A)\times\cdots=0$ too.