Grain growth is a well-known and complex phenomenon occurring during annealing of all polycrystalline materials. Its numerical modeling is a complex task when anisotropy sources such as grain orientation and grain boundary inclination have to be taken into account. This article presents the application of the front-tracking methodology ToRealMotion introduced in previous works, to the context of anisotropic grain boundary motion at the mesoscopic scale. The new formulation of boundary migration can take into account any source of anisotropy both at grain boundaries as well as at multiple junctions (MJs) (intersection point of three or more grain boundaries). Special attention is given to the decomposition of high-order MJs for which an algorithm is proposed based on local grain boundary energy minimisation. Numerical tests are provided using highly heterogeneous configurations, and comparisons with a recently developed Finite-Element Level-Set (FE-LS) approach are given. Finally, the computational performance of the model will be studied comparing the CPU-times obtained with the same model but in an isotropic context.