$U$ bisected by all hyperplanes $\implies$ $U$ symmetric? Let $U \subset \mathbb{R}^n$ be a bounded open convex set, such that every hyperplane passing through the origin divides $U$ into two sets of equal volume ($n$-dimensional Lebesgue measure?). 

Does it follow that $U$ is symmetric with respect to the origin, i.e. $U = -U$?

I'm fairly sure this holds for $n=2$: If two distinct hyperplanes (lines) through the origin divide $U$ into four sets $A, B, C, D$, then we have 
\begin{align*}
\text{vol}(U)/2 &= \text{vol}(A)+\text{vol}(B) = \text{vol}(C)+\text{vol}(D) \\
                &= \text{vol}(B)+\text{vol}(C) = \text{vol}(D)+\text{vol}(A)
\end{align*}
Hence $\text{vol}(A) = \text{vol}(C)$, $\text{vol}(B) = \text{vol}(D)$. The boundary of $U$ is a simple closed curve, and a limiting argument using this equality of volumes implies that opposite points on this curve ("opposite" meaning the line between them passes through the origin) have the same distance from the origin. Therefore the boundary of $U$ is symmetric, so $U$ is symmetric.
For $n > 2$, if $U$ is divided by $n$ hyperplanes, then opposite volumes are not necessarily equal, so the above approach does not apply.
 A: (Leaving as CW since it does not resolve the question, but just indicates some possible points of view.)
Walsh functions
The fact that the argument that the OP gave works in 2 but higher dimensions is related to properties of the Walsh functions. In $n$ dimensions with $n$ hyperplanes we divide the body into $2^n$ regions. The regions can be numbered by binary numbers as follows: the $k$th bit of the number is $1$ if it the region is "above" the $k$th hyperplane, and $0$ if the region is "below" the $k$th hyperplane. 
Now let $V$ be the vector of length $2^n$ whose entries are the volumes of the corresponding regions, ordered according to the numbering above. The reflection condition gives $n$ equations of the form 
$$ A_k \cdot V = 0 \tag{W1}$$
where $A_k$ is a vector of length $2^n$, whose entries take value $\pm1$ corresponding to whether the corresponding region lies "above" or "below" the $k$th hyperplane. 
It is clear that each $A_k$ is one of the Walsh functions for $2^n$ dimension. 
The method that works in 2 dimensions as outlined in the OP is the claim that the equation (W1) implies
$$ (1,0\ldots, 0, -1) \cdot V = 0 \tag{W2}$$
And this is true in 2 dimension since of the four Walsh functions, namely
$$ (1,1,1,1) \\ (1,1,-1,-1) \\ (1,-1,-1,1) \\ (1,-1,1,-1) $$
the first and third are orthogonal to $(1,0,0,-1)$ and the second and fourth are exactly the ones corresponding to $A_1$ and $A_2$ in (W1). Using that the Walsh functions form an orthogonal bases for $\mathbb{R}^{2^n}$ we see that this allows us to reconstruct $(1,0,0,-1)$ as a superposition of $A_1$ and $A_2$. 
In $>2$ dimensions this is no longer true, simply from a counting argument. $(1,0,0,\ldots,0,-1)$ is orthogonal to exactly half of the Walsh functions, so you need knowledge of exactly half of the Walsh functions to reconstruct it. However, for $n > 2$ we have $2^{n-1} > n$ so (W2) cannot be derived as a consequence of (W1). 
A related problem
Instead of considering the volume of a convex body, we can consider the following set-up: 

Problem 2 Let $f:\mathbb{S}^{n-1} \to \mathbb{R}$ be a Lebesgue-integrable function defined on the sphere. Given a hyperplane $\Pi$ through the origin, we denote by $H_\pm$ the two hemispheres. Suppose that for every $\Pi$ it is true that $\int_{H_+} f \mathrm{d}\sigma = \int_{H_-} f \mathrm{d}\sigma$, is it true that $f$ is an even function? ($f(-x) = f(x)$)

A positive answer to the above would give a positive answer to the question in the OP, since the volume for the convex body can be computed by performing the radial integral first and then integrating over the sphere, which reduces it to a question like  the above. 
Of course, for Lebesgue integrable functions the claim is not true, at least not without replacing "$f$ is even" by "$f$ is even almost everywhere." 
We note then the proof given in the OP also solves the $n = 2$ case of Problem 2 if we assume $f$ is continuous; this extends the result from "convex" bodies to "star-shaped bodies" such that the distance from the boundary to the center is a continuous function on the sphere. 
Harmonic analysis
One can replace the Walsh functions approach by a harmonic analysis approach, considering Problem 2. Suppose $f$ is $L^2$ integrable on the sphere, we know that it can be decomposed into spherical harmonics. Using the symmetry properties of the spherical harmonics, Problem 2 would be resolved if we can show that the spherical harmonic coefficients for $f$ vanishes whenever angular momentum parameter $\ell$ is odd. 
One method to approach this is to show that each spherical harmonic can be reconstructed from the hemispherical integrals by superposition. More precisely, let $v\in \mathbb{S}^{n-1}$ and let $\chi_v(x) = \mathrm{sgn}(x\cdot v)$, this function is positive on one hemisphere and negative on the other. Our goal then is to reconstruct the spherical harmonics by finding for each harmonic $Y$ a corresponding $\rho$ such that
$$ Y(x) = \int \rho(v) \chi_v(x) \mathrm{d}v \tag{H1}$$
This is the continuous version of what happened in the Walsh functions section. 
It turns out that in 2 dimensions solving (H1) is pretty easy: to reconstruct the harmonic $\sin(k \theta)$ for $\theta\in [0,2\pi)$ we can let $\rho$ be $\frac1{4k} \cos (k\theta)$. But as to whether this can be done in higher dimensions: I have no idea. 
