This paper discusses quaternion $L^p$ geometric weighting averaging working on the multiplicative Lie group of nonzero quaternions $\mathbb{H}^{*}$, endowed with its natural bi-invariant Riemannian metric. Algorithms for computing the Riemannian $L^p$ center of mass of a set of points, with $1 \leq p \leq \infty$ (i.e., median, mean, $L^p$ barycenter and minimax center), are particularized to the case of $\mathbb{H}^{*}$.Two different approaches are considered. The first formulation is based on computing the logarithm of quaternions which maps them to the Euclidean tangent space at the identity $\mathbf{1}$, associated to the Lie algebra of $\mathbb{H}^{*}$. In the tangent space, Euclidean algorithms for $L^p$ center of mass can be naturally applied. The second formulation is a family of methods based on gradient descent algorithms aiming at minimizing the sum of quaternion geodesic distances raised to power $p$. These algorithms converges to the quaternion Fr\'{e}chet-Karcher barycenter ($p=2$), the quaternion Fermat-Weber point ($p=1$) and the quaternion Riemannian 1-center ($p=+\infty$). Besides giving explicit forms of these algorithms, their application for quaternion image processing is shown by introducing the notion of quaternion bilateral filtering