When implementing the vanilla algorithm of Vicsek model from Vicsek’s 1995 PRL paper in 3D, one need to rotate a lot of vectors. This post introduced my idea of generating such noise.
2D and 3D
In the original paper, the noise is implemented by rotating the targeted direction by a random direction. I find the math equation in this paper being very good addressing such process, as the following.
\[\mathbf{s}_{i}^{t+\Delta t}=\mathcal{R}_{\eta}\left[\frac{\sum_{j} n_{i j}^{t} \mathbf{s}_{j}^{t}}{\left\|\sum_{j} n_{i j}^{t} \mathbf{s}_{j}^{t}\right\|}\right]\]The operator $\mathcal{R}_\eta(\mathbf{s})$ is responsible for rotating vector $\mathbf{s}$ around for a magnitude of $\eta$ randomly.
In the 2D case, if the angle of vector $\mathbf{s}$ is $\theta$, then the operator $\mathcal{R}_\eta(\mathbf{s})$ actually does this.
\[\theta \rightarrow \theta + \text{U}(-\pi \eta, \pi\eta)\]where $U$ corresponds to a uniform distribution.
In 3D however, things are getting tricker. The following graph shows the effect of $\mathcal{R}_\eta$. The desired effect is to rotate angle $s$ randomly inside the pink region.
The next chapter will discuss the way to perform this operation.
Generating 3D noise
The practical way to generate scalar noise in 3D is this,
- Generate uniformly distributed noise as if vector $\mathbf{s}$ is (0, 0, 1).
- Rotate the generated noise to the direction of $\mathbf{s}$
Random Vectors
The aim is to generate random vectors pointing uniformly to the pink cap.
The python code that generating such noise is this
noise_phi = np.random.uniform(-np.pi, np.pi, N)
noise_z = np.random.uniform(1 - 2 * eta, 1, N)
noise_r = np.sqrt(1 - noise_z**2)
noise_x = noise_r * np.cos(noise_phi)
noise_y = noise_r * np.sin(noise_phi)
noise = np.array((noise_x, noise_y, noise_z)).T
The not very obvious code is this line.
noise_z = np.random.uniform(1 - 2 * eta, 1, N)
Here, $\eta$ (eta) is a measure of noise level ranging from $0$ to $1$. When $\eta = 1$ the pink cap would be the entire unit sphere, and the cap will vanish if $\eta = 0$.
Embarrassingly I learned how to pick up points on a sphere just recently at the age of 28, with this extremely helpful stackexchange answer. Basically the idea is if I need the uniform distribution on the pink cap, I will need two uniform distributions.
\[\begin{aligned} \theta &\sim \text{U}(-\pi, \pi)\\ z &\sim \text{U}(1-h, 1)\\ &\sim \text{U}(1-2\eta, 1) \end{aligned}\]The relationship between $\eta$ (eta) and $h$, $h \sim 2\eta$, comes from the area of the pink cap inside a unit sphere.
where $h$ takes the range between $0$ and $2$.
Rotation
Now we have noise vectors around $\mathbf{s}’$, the next step is to find a rotation matrix $\mathbf{R}$ that performs the following transformation.
\[\bf{s = R \cdot s'}\]
following this stackexchange answer, it is intuitive to break the rotation into the following parts.
- Transform the basis from $\bf{(x, y, z)}$ to $\bf{(u, v, w)}$.
- Perform a simple rotation.
- Transform back from $\bf{(u, v, w)}$ to original basis $\bf{(x, y, z)}$.
The key observation to calculate this rotation is that such rotation is very straightforward in the following $\bf{(u, v, w)}$ basis. It is merely a rotation along $\bf{w}$ axis.
- $\bf{u}$: the vector projection of $\bf{s}$ on $\bf{s}’$
- $\bf{v}$: the vector rejection of $\bf{s}$ on $\bf{s}’$
- $\bf{w}$: the cross product of $\bf{s}$ and $\bf{s}’$
The rotation angle $\theta$ in the $\bf(u,v,w)$ basis is calculated by $\arccos(\bf{s \cdot s’})$, and the rotation matrix is simply this.
\[\begin{pmatrix} \cos\theta & -\sin\theta & 0\\ \sin\theta & \cos\theta & 0\\ 0 & 0 & 1 \end{pmatrix}\]A step-by-step operation tutorial is in this stackexchange answer.