I've long been intrigued by this surface which tightly hugs each axis, extending to infinity but with finite volume. But integrating this formula is beyond my powers. Any suggestions on how to do this, or demonstrations of the full answer? I suspect spherical coordinates might help, but that seems to make the formulas even more complex and beyond my capacity.
An upper bound on volume: Consider the surface near the x-axis only, for x > 1. The cross-section around the axis should be within the circle $y^2 + z^2 = 1/x^2$ [Gabriel's Horn, as David K helpfully pointed out below], and yet near it, since the dropped term, $y^2 z^2 / x^2$, is << $1/x^2$.
This circle has area $\pi/x^4$; integrating from 1 to infinity gives an upper bound of volume around all axes of $6\pi$ plus the area within 1 unit of the origin (max. 8 there, of course), for a total <27. The graph of the cross-section still looks nearly circular at x=1, so I'm expecting a more exact answer might be in the 24-26 range.
Slight progress: the smallest cube enclosed by the surface has vertices at x=y=z = $1/\sqrt[4]3$, so its volume is $8/\sqrt[4]{27}$. The six arms fully separate at axis values >1. In between these values one would have to integrate a complex cross section of squares with rounded corners--the square sides separating the parts closest to each axis to get a partial volume, then multiplying by 6 just as for each arm segment to get a total volume. But both this and the infinite part of the arm requires the anti-differential of $\sqrt{((1-x^2y^2)/(x^2+y^2))}$.
Latest update: I found that wolfram alpha has a double integral calculator widget; it doesn't show an anti-derivative for the above function of z, but presumably did a numerical integration and delivered the volume for the curve (x: 0-infinity; y: 0-1/x) of 3.24099; multiplying by 8 gives 25.928, in the ballpark of my earlier estimate. It would still be interesting to know if it has some more exact formulation (even if this involves square roots or other complex terms).
Take $\mathbb{E}$ to be the region of space given by $$\mathbb{E}=\{(x,y,z)\in\mathbb(0,\infty)^3|(xy)^2+(xz)^2+(yz)^2 < 1\}.$$ If we induce the change of coordinates $(u,v,w)=(xy,xz,yz)$ we can say $$\int_{\mathbb{E}}dV=\int_{\{u^2+v^2+w^2 < 1\}\cap(0,\infty)^3}\left| \frac{\partial(x,y,z)}{\partial(u,v,w)}\right|dV=\int_0^{1}\int_0^{\sqrt{1-u^2}}\int_0^{\sqrt{1-u^2-v^2}}\Bigg(\frac{1}{2\sqrt{uvw}}\Bigg)dwdvdu$$ After converting to cylindrical coordinates and integrating we get $$\int_{\mathbb{E}}dV=\Bigg(\int_0^{\pi/2}\frac{dx}{\sqrt{\sin(x)\cos(x)}}\Bigg)\cdot \Bigg(\int_0^1\sqrt[4]{1-x^2}dx\Bigg)$$ Hence we can express the volume of our solid as $$8\cdot\Bigg(\int_0^{\pi/2}\frac{dx}{\sqrt{\sin(x)\cos(x)}}\Bigg)\cdot \Bigg(\int_0^1\sqrt[4]{1-x^2}dx\Bigg)$$ which reduces the problem to calculating two real-valued integrals.