I am working with voxel volumes which are addressed using spherical coordinates. That is, I operate on voxel elements $V(r, \theta, \phi)$ where $r$ is range, $\theta$ is the azimuth, and $\phi$ is the inclination. I am however running into issues with the fact that the angular ranges wrap around.
To make the problem a bit easier to follow, let's look at only the azimuth. The azimuth $\theta$ is defined to lie on the interval $\theta \in [-\pi, \pi]$ and wraps around such that $\theta + 2\pi = \theta$. I now want to map any azimuth value to a voxel index, which is the same as gridding the data. My current implementation is to simply map this as
$$ i(\theta) = v_n \cdot \left\lfloor \frac{\theta - \theta_0}{v_s} \right\rfloor $$
where $v_n$ is the number of voxels (grid points), $\theta_0$ is the origin such that $i(\theta_0) = 0$, and $v_s$ is the voxel scale (grid size).
This works, but gives a lot of empty voxels in case the data is defined near $\pm \pi$. For example, if my data range is $\theta \in [-\pi, -\pi + a] \cup [\pi-b, \pi]$ (the same as $[\pi - a, \pi + b]$ before wrapping) then the angular range is $a + b$, but I would have to have voxels created for the entire $2 \pi$ range.
So my question is how do I best define the mapping between an azimuth angle and its corresponding voxel index?, with the constraints that
- I use as few voxels as possible for the data
- The mapping is invertible so that I can get the spherical coordinate of the voxel center, given its index ($\theta_c(i) = \ldots$)
My only attempt so far is to detect whether the data is best represented by an azimuth range of $[-\pi, \pi]$ or $[0, 2\pi]$, and use two different mappings. But this seems error prone and inelegant.