diff --git a/NEWS.md b/NEWS.md
index 6d5d6835d..821577357 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -43,6 +43,9 @@
   The new features are more extensively described in help pages
   `?plot.cca`, `?ordiplot` and `?biplot.rda`.
 
+* `rda` and `cca` return centroids for factor levels also when
+  called without formula, for instance `cca(dune, dune.env)`.
+
 ## Bug Fixes
 
 * `summary.ordihull` failed if input data were not two-dimensional.
diff --git a/R/capscale.R b/R/capscale.R
index 963202318..ab4874f56 100644
--- a/R/capscale.R
+++ b/R/capscale.R
@@ -134,17 +134,9 @@
             sol$CCA$v[] <- NA
         sol$colsum <- NA
     }
-    if (!is.null(sol$CCA) && sol$CCA$rank > 0)
-        sol$CCA$centroids <- centroids.cca(sol$CCA$wa, d$modelframe)
-    if (!is.null(sol$CCA$alias))
-        sol$CCA$centroids <- unique(sol$CCA$centroids)
-    if (!is.null(sol$CCA$centroids)) {
-        rs <- rowSums(sol$CCA$centroids^2)
-        sol$CCA$centroids <- sol$CCA$centroids[rs > 1e-04, ,
-                                               drop = FALSE]
-        if (nrow(sol$CCA$centroids) == 0)
-            sol$CCA$centroids <- NULL
-    }
+    ## centroids
+    sol$CCA$centroids <- getCentroids(sol, d$modelframe)
+
     sol$call <- match.call()
     sol$terms <- terms(formula, "Condition", data = data)
     sol$terminfo <- ordiTerminfo(d, data)
diff --git a/R/cca.default.R b/R/cca.default.R
index 482e86741..63794d810 100644
--- a/R/cca.default.R
+++ b/R/cca.default.R
@@ -8,8 +8,10 @@
         stop("function cannot be used with (dis)similarities")
     X <- as.matrix(X)
     if (!is.null(Y)) {
-        if (is.data.frame(Y) || is.factor(Y))
+        if (is.data.frame(Y) || is.factor(Y)) { # save Y for centroids
+            mframe <- as.data.frame(Y) # can be a single factor
             Y <- model.matrix(~ ., as.data.frame(Y))[,-1,drop=FALSE]
+        }
         Y <- as.matrix(Y)
     }
     if (!is.null(Z)) {
@@ -27,6 +29,10 @@
         X <- X[, !tmp, drop = FALSE]
     }
     sol <- ordConstrained(X, Y, Z, method = "cca")
+    ## mframe exists only if function was called as cca(X, mframe)
+    if (exists("mframe"))
+        sol$CCA$centroids <- getCentroids(sol, mframe)
+
     if (exists("exclude.spec")) {
         if (!is.null(sol$CCA$v))
             attr(sol$CCA$v, "na.action") <- exclude.spec
diff --git a/R/cca.formula.R b/R/cca.formula.R
index 9ca7b6c9b..648faca99 100644
--- a/R/cca.formula.R
+++ b/R/cca.formula.R
@@ -10,21 +10,7 @@
     d <- ordiParseFormula(formula, data = data, na.action = na.action,
                           subset = substitute(subset))
     sol <- cca.default(d$X, d$Y, d$Z)
-    if (!is.null(sol$CCA) && sol$CCA$rank > 0) {
-        centroids <- centroids.cca(sol$CCA$wa, d$modelframe,
-                                   sol$rowsum)
-        if (!is.null(sol$CCA$alias))
-            centroids <- unique(centroids)
-        ## See that there really are centroids
-        if (!is.null(centroids)) {
-            rs <- rowSums(centroids^2)
-            centroids <- centroids[rs > 1e-04,, drop = FALSE]
-            if (length(centroids) == 0)
-                centroids <- NULL
-        }
-        if (!is.null(centroids))
-            sol$CCA$centroids <- centroids
-    }
+    sol$CCA$centroids <- getCentroids(sol, d$modelframe)
     ## replace cca.default call
     call <- match.call()
     call[[1]] <- as.name("cca")
diff --git a/R/centroids.cca.R b/R/centroids.cca.R
index 29ba2b079..4b1626f3b 100644
--- a/R/centroids.cca.R
+++ b/R/centroids.cca.R
@@ -9,7 +9,7 @@
     mf <- mf[, facts, drop = FALSE]
     ## Explicitly exclude NA as a level
     mf <- droplevels(mf, exclude = NA)
-    if (missing(wt))
+    if (missing(wt) || is.null(wt))
         wt <- rep(1, nrow(mf))
     ind <- seq_len(nrow(mf))
     workhorse <- function(x, wt)
@@ -36,3 +36,30 @@
     colnames(out) <- colnames(x)
     out
 }
+
+### cca.centroids is used by all constrained ordination methods and
+### factorfit (via envfit). All constrained ordinations sanitize the
+### results is the same way, and instead of having the same code in
+### all functions, let's have it here.
+
+## @param ord ordConstrained result object.
+## @param mframe Data frame, possibly with factors for which centroids
+## are needed.
+
+`getCentroids` <-
+    function(ord, mframe)
+{
+    if (is.null(ord$CCA) || ord$CCA$rank < 1)
+        return(NULL)
+    centroids <- centroids.cca(ord$CCA$u, mframe, ord$rowsum)
+    if (!is.null(ord$CCA$alias))
+        centroids <- unique(centroids)
+    ## See that there really are centroids
+    if (!is.null(centroids)) {
+        rs <- rowSums(centroids^2)
+        centroids <- centroids[rs > 1e-04,, drop = FALSE]
+        if (length(centroids) == 0)
+            centroids <- NULL
+    }
+    centroids
+}
diff --git a/R/dbrda.R b/R/dbrda.R
index 1af558ce8..7545fd4f4 100644
--- a/R/dbrda.R
+++ b/R/dbrda.R
@@ -115,18 +115,8 @@
                                        drop = FALSE]
         sol$CA$u <- sol$CA$u[, seq_len(sol$CA$poseig), drop = FALSE]
     }
-    if (!is.null(sol$CCA) && sol$CCA$rank > 0)
-        sol$CCA$centroids <-
-            centroids.cca(sol$CCA$u, d$modelframe)
-    if (!is.null(sol$CCA$alias))
-        sol$CCA$centroids <- unique(sol$CCA$centroids)
-    if (!is.null(sol$CCA$centroids)) {
-        rs <- rowSums(sol$CCA$centroids^2)
-        sol$CCA$centroids <- sol$CCA$centroids[rs > 1e-04, ,
-                                               drop = FALSE]
-        if (nrow(sol$CCA$centroids) == 0)
-            sol$CCA$centroids <- NULL
-    }
+    sol$CCA$centroids <- getCentroids(sol, d$modelframe)
+
     sol$call <- match.call()
     sol$terms <- terms(formula, "Condition", data = data)
     sol$terminfo <- ordiTerminfo(d, data)
diff --git a/R/rda.default.R b/R/rda.default.R
index 4dd29134d..1f17f6f3a 100644
--- a/R/rda.default.R
+++ b/R/rda.default.R
@@ -8,8 +8,10 @@
         stop("function cannot be used with (dis)similarities")
     X <- as.matrix(X)
     if (!is.null(Y)) {
-        if (is.data.frame(Y) || is.factor(Y))
+        if (is.data.frame(Y) || is.factor(Y)) { # save Y for centroids
+            mframe <- as.data.frame(Y) # can be a single factor
             Y <- model.matrix(~ ., as.data.frame(Y))[,-1,drop=FALSE]
+        }
         Y <- as.matrix(Y)
     }
     if (!is.null(Z)) {
@@ -19,6 +21,9 @@
     }
 
     sol <- ordConstrained(X, Y, Z, arg = scale, method = "rda")
+    ## mframe exists only if function was called rda(X, mframe)
+    if (exists("mframe"))
+        sol$CCA$centroids <- getCentroids(sol, mframe)
 
     call <- match.call()
     call[[1]] <- as.name("rda")
diff --git a/R/rda.formula.R b/R/rda.formula.R
index 439781753..aec9db9ff 100644
--- a/R/rda.formula.R
+++ b/R/rda.formula.R
@@ -11,19 +11,7 @@
     d <- ordiParseFormula(formula, data = data, na.action = na.action,
                           subset = substitute(subset))
     sol <- rda.default(d$X, d$Y, d$Z, scale)
-    if (!is.null(sol$CCA) && sol$CCA$rank > 0) {
-        centroids <- centroids.cca(sol$CCA$wa, d$modelframe)
-        if (!is.null(sol$CCA$alias))
-            centroids <- unique(centroids)
-        if (!is.null(centroids)) {
-            rs <- rowSums(centroids^2)
-            centroids <- centroids[rs > 1e-04,, drop = FALSE]
-            if (length(centroids) == 0)
-                centroids <- NULL
-        }
-        if (!is.null(centroids))
-            sol$CCA$centroids <- centroids
-    }
+    sol$CCA$centroids <- getCentroids(sol, d$modelframe)
     ## replace rda.default call
     call <- match.call()
     call[[1]] <- as.name("rda")
diff --git a/man/plot.cca.Rd b/man/plot.cca.Rd
index 5346e6a72..8c50cc932 100644
--- a/man/plot.cca.Rd
+++ b/man/plot.cca.Rd
@@ -173,9 +173,7 @@
   function uses a simple heuristics for adjusting the unit-length arrows
   to the current plot area, but the user can give the expansion factor
   in \code{arrow.mul}.  With \code{display="cn"} the centroids of levels
-  of \code{\link{factor}} variables are displayed (these are available
-  only if there were factors and a formula interface was used in
-  \code{\link{cca}} or \code{\link{rda}}). With this option continuous
+  of \code{\link{factor}} variables are displayed. With this option continuous
   variables still are presented as arrows and ordered factors as arrows
   and centroids. With \code{display = "reg"} arrows will be drawn for
   regression coefficients (a.k.a. canonical coefficients) of constraints