diff --git a/DESCRIPTION b/DESCRIPTION
index 4e01f4f..76ef035 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -32,7 +32,7 @@ Imports:
     ggrepel,
     ggforce,
     tidyr
-RoxygenNote: 7.2.0
+RoxygenNote: 7.2.3
 Depends: 
     R (>= 2.10)
 Suggests: 
diff --git a/R/bind.R b/R/bind.R
index b07608d..0c59f28 100644
--- a/R/bind.R
+++ b/R/bind.R
@@ -50,23 +50,15 @@ bind_theoretical <- function(dt, matrix, index, raw_data) {
 #' @return a tibble object containing both the searched and random bases
 #' @export
 bind_random <- function(dt, n = 500, seed = 1) {
-  p <- nrow(dt$basis[[1]]) * ncol(dt$basis[[1]])
-  ncol <- nrow(dt$basis[[1]])
-
-  fix_matrix <- function(dt) {
-    dt <- t(as.matrix(dt, nrow = p, ncol = ncol))
-    rownames(dt) <- NULL
-    dt
-  }
+  m <- dt$basis[[1]]
 
   set.seed(seed)
-  n_geozoo <- p * n
-  suppressWarnings(sphere_basis <- geozoo::sphere.hollow(p, n_geozoo)$points %>%
-    tibble::as_tibble() %>%
-    dplyr::nest_by(id = dplyr::row_number()) %>%
-    dplyr::ungroup() %>%
-    dplyr::mutate(basis = purrr::map(.data$data, fix_matrix)) %>%
-    dplyr::select(.data$basis))
+  sphere_basis <- purrr::map_dfr(
+    1:n, ~tibble::tibble(
+      id = .x,
+      basis = list(tourr::basis_random(n = nrow(m), d = ncol(m)))
+      )
+    )
 
   sphere_points <- sphere_basis %>%
     dplyr::mutate(
@@ -86,6 +78,7 @@ bind_random <- function(dt, n = 500, seed = 1) {
 #'
 #' @param basis a matrix returned by \code{get_basis_matrix()}
 #' @param n numeric; the number of random bases to generate in each dimension by geozoo
+#' @param d numeric; dimension of the basis, d = 1, 2, ...
 #' @param front logical; if the random bases should be bound before or after the original bases
 #' @param seed numeric; a seed for generating reproducible random bases from geozoo
 #' @examples
@@ -95,18 +88,18 @@ bind_random <- function(dt, n = 500, seed = 1) {
 #' @family bind
 #' @return a matrix containing both the searched and random bases
 #' @export
-bind_random_matrix <- function(basis, n = 500, front = FALSE, seed = 1) {
-  p <- ncol(basis)
-  n_geozoo <- p * n
-  set.seed(seed)
-  sphere_basis <- geozoo::sphere.hollow(p, n_geozoo)$points
-  colnames(sphere_basis) <- colnames(basis)
+bind_random_matrix <- function(basis, n = 500, d = 1, front = FALSE, seed = 1) {
 
+  r <- ncol(basis)/d
+  set.seed(seed)
+  random_basis <- do.call(
+    rbind, purrr::map(1:n, ~t(tourr::basis_random(n = r, d = d)))
+    )
 
   if (front) {
-    out <- sphere_basis %>% rbind(basis)
+    out <- random_basis %>% rbind(basis)
   } else {
-    out <- basis %>% rbind(sphere_basis)
+    out <- basis %>% rbind(random_basis)
   }
 
   return(out)
diff --git a/R/explore_space.R b/R/explore_space.R
index 5c5b000..d1e22f4 100644
--- a/R/explore_space.R
+++ b/R/explore_space.R
@@ -239,13 +239,13 @@ compute_pca <- function(dt, group = NULL, random = TRUE, flip = TRUE, ...) {
 #'   \item{\code{prep_space_tour()}}{a list containing various components needed for producing the animation}
 #' }
 #' @export
-explore_space_tour <- function(...) {
+explore_space_tour <- function(..., axes = "bottomleft") {
   prep <- prep_space_tour(...)
 
   tourr::animate_xy(prep$basis,
                     col = prep$col, cex = prep$cex, pch = prep$pch,
                     edges = prep$edges, edges.col = prep$edges_col,
-                    axes = "bottomleft"
+                    axes = axes
   )
 }
 
@@ -253,14 +253,17 @@ explore_space_tour <- function(...) {
 #' @param dt a data object collected by the projection pursuit guided tour optimisation in \code{tourr}
 #' @param group the variable to label different runs of the optimiser(s)
 #' @param flip logical; if the sign flipping need to be performed
+#' @param n_random numeric; the number of random basis to generate
 #' @param color the variable to be coloured by
 #' @param rand_size numeric; the size of random points
+#' @param rand_color character; the color hex code for random points
 #' @param point_size numeric; the size of points searched by the optimiser(s)
 #' @param end_size numeric; the size of end points
 #' @param theo_size numeric; the size of theoretical point(s)
 #' @param theo_shape numeric; the shape symbol in the basic plot
 #' @param theo_color character; the color of theoretical point(s)
 #' @param palette the colour palette to be used
+#' @param axes see [tourr::animate_xy()]
 #' @param ... other argument passed to \code{tourr::animate_xy()} and \code{prep_space_tour()}
 #' @examples
 #' explore_space_tour(dplyr::bind_rows(holes_1d_better, holes_1d_geo),
@@ -268,9 +271,10 @@ explore_space_tour <- function(...) {
 #' )
 #' @rdname explore_space_tour
 #' @export
-prep_space_tour <- function(dt, group = NULL, flip = FALSE,
-                            color = NULL, rand_size = 1, point_size = 1.5, end_size = 5,
-                            theo_size = 3, theo_shape = 17, theo_color = "black",
+prep_space_tour <- function(dt, group = NULL, flip = FALSE, n_random = 2000,
+                            color = NULL, rand_size = 1, rand_color = "#D3D3D3",
+                            point_size = 1.5, end_size = 5, theo_size = 3,
+                            theo_shape = 17, theo_color = "black",
                             palette = botanical_palettes$fern, ...) {
   if (rlang::quo_is_null(dplyr::enquo(color))) {
     message("map method to color")
@@ -287,7 +291,9 @@ prep_space_tour <- function(dt, group = NULL, flip = FALSE,
     basis <- flip$basis %>% bind_random_matrix(front = TRUE)
   } else{
     flip = list(dt = dt)
-    basis <- dt %>% get_basis_matrix() %>% bind_random_matrix(front = TRUE)
+    basis <- dt %>%
+      get_basis_matrix() %>%
+      bind_random_matrix(n = n_random, front = TRUE)
   }
 
   n_rand <- nrow(basis) - nrow(dt)
@@ -309,7 +315,7 @@ prep_space_tour <- function(dt, group = NULL, flip = FALSE,
   edges_col <- palette[as.factor(edges_dt %>% dplyr::pull({{ color }}))]
 
   col <- c(
-    rep("#D3D3D3", n_rand),
+    rep(rand_color, n_rand),
     palette[as.factor(dt %>% dplyr::pull({{ color }}))]
   )
   cex <- c(
diff --git a/man/bind_random_matrix.Rd b/man/bind_random_matrix.Rd
index 77cd8f7..5982efd 100644
--- a/man/bind_random_matrix.Rd
+++ b/man/bind_random_matrix.Rd
@@ -4,13 +4,15 @@
 \alias{bind_random_matrix}
 \title{Bind random bases in the projection bases space as a matrix}
 \usage{
-bind_random_matrix(basis, n = 500, front = FALSE, seed = 1)
+bind_random_matrix(basis, n = 500, d = 1, front = FALSE, seed = 1)
 }
 \arguments{
 \item{basis}{a matrix returned by \code{get_basis_matrix()}}
 
 \item{n}{numeric; the number of random bases to generate in each dimension by geozoo}
 
+\item{d}{numeric; dimension of the basis, d = 1, 2, ...}
+
 \item{front}{logical; if the random bases should be bound before or after the original bases}
 
 \item{seed}{numeric; a seed for generating reproducible random bases from geozoo}
diff --git a/man/explore_space_tour.Rd b/man/explore_space_tour.Rd
index 9e1901d..a23b602 100644
--- a/man/explore_space_tour.Rd
+++ b/man/explore_space_tour.Rd
@@ -5,14 +5,16 @@
 \alias{prep_space_tour}
 \title{Plot the grand tour animation of the bases space in high dimension}
 \usage{
-explore_space_tour(...)
+explore_space_tour(..., axes = "bottomleft")
 
 prep_space_tour(
   dt,
   group = NULL,
   flip = FALSE,
+  n_random = 2000,
   color = NULL,
   rand_size = 1,
+  rand_color = "#D3D3D3",
   point_size = 1.5,
   end_size = 5,
   theo_size = 3,
@@ -25,16 +27,22 @@ prep_space_tour(
 \arguments{
 \item{...}{other argument passed to \code{tourr::animate_xy()} and \code{prep_space_tour()}}
 
+\item{axes}{see [tourr::animate_xy()]}
+
 \item{dt}{a data object collected by the projection pursuit guided tour optimisation in \code{tourr}}
 
 \item{group}{the variable to label different runs of the optimiser(s)}
 
 \item{flip}{logical; if the sign flipping need to be performed}
 
+\item{n_random}{numeric; the number of random basis to generate}
+
 \item{color}{the variable to be coloured by}
 
 \item{rand_size}{numeric; the size of random points}
 
+\item{rand_color}{character; the color hex code for random points}
+
 \item{point_size}{numeric; the size of points searched by the optimiser(s)}
 
 \item{end_size}{numeric; the size of end points}
diff --git a/tests/testthat/test-bind.R b/tests/testthat/test-bind.R
index 637e404..9d4bbe6 100644
--- a/tests/testthat/test-bind.R
+++ b/tests/testthat/test-bind.R
@@ -40,5 +40,5 @@ data <- get_basis_matrix(holes_1d_geo)
 
 test_that("bind_random_matrix", {
   expect_true(bind_random_matrix(data) %>% is.matrix())
-  expect_equal(bind_random_matrix(data) %>% nrow(), nrow(data) + ncol(data) * 500)
+  expect_equal(bind_random_matrix(data) %>% nrow(), nrow(data) + 500)
 })