diff --git a/DESCRIPTION b/DESCRIPTION
index 2c1a82f..92bc807 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,5 +35,5 @@ Suggests:
 License: AGPL-3
 LazyData: true
 ByteCompile: true
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 VignetteBuilder: knitr
diff --git a/NEWS.md b/NEWS.md
index 76cccc4..d60873c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,7 @@
 # foto 1.2
 
 * fix plotting routine
+* flexible constraints on harmonics
 * adding new vignette - theoretical basis
 
 # foto 1.1
diff --git a/R/foto.R b/R/foto.R
index d123358..4f995c2 100644
--- a/R/foto.R
+++ b/R/foto.R
@@ -7,6 +7,7 @@
 #' @param x an image file, or single or multi-layer SpatRaster
 #' (RGB or otherwise), multi-layer data are averaged to a single layer
 #' @param window_size a moving window size in pixels (default = 61 pixels)
+#' @param harmonics number of harmonics to consider (29 by default)
 #' @param method zones (for discrete zones) or mw for a moving window
 #' approach
 #' @param pca execute PCA, \code{TRUE} or \code{FALSE}. If \code{FALSE} only
@@ -41,13 +42,16 @@
 #' print(names(output))
 #' }
 #'
-foto <- function(x,
-                 window_size = 61,
-                 method = "zones",
-                 norm_spec = FALSE,
-                 high_pass = TRUE,
-                 pca = TRUE,
-                 plot = FALSE) {
+foto <- function(
+    x,
+    window_size = 61,
+    harmonics = 29,
+    method = "zones",
+    norm_spec = FALSE,
+    high_pass = TRUE,
+    pca = TRUE,
+    plot = FALSE
+    ) {
   
   # get the current enviroment
   env <- environment()
@@ -99,11 +103,7 @@ foto <- function(x,
 
   # use assign
   # define output matrix and global increment, i
-  if (window_size / 2 < 29) {
-    output <- matrix(0, cells + 1, window_size / 2)
-  } else {
-    output <- matrix(0, cells, 29)
-  }
+  output <- matrix(0, cells, harmonics)
 
   # see rspectrum function above
   i <- 0
@@ -117,6 +117,7 @@ foto <- function(x,
           rspectrum(
             x = x,
             w = window_size,
+            hm = harmonics,
             n = norm_spec,
             h = high_pass,
             env = env
@@ -138,6 +139,7 @@ foto <- function(x,
           rspectrum(
             x = x,
             w = window_size,
+            hm = harmonics,
             n = norm_spec,
             env = env
           )
diff --git a/R/rspectrum.R b/R/rspectrum.R
index 5e121ba..4fb37e8 100644
--- a/R/rspectrum.R
+++ b/R/rspectrum.R
@@ -4,6 +4,7 @@
 #'
 #' @param x a square matrix
 #' @param w a moving window size
+#' @param hm harmonics to consider as features
 #' @param n normalize, bolean \code{TRUE} or \code{FALSE}
 #' @param h high pass filter on the two first spectra values
 #'  set to 0, limits the influence of low frequency components
@@ -13,12 +14,16 @@
 #' @return Returns a radial spectrum values for the image used
 #' in order to classify texture using a PCA (or other) analysis.
 
-rspectrum <- function(x,
-                      w,
-                      n = TRUE,
-                      h = TRUE,
-                      env,
-                      ...) {
+rspectrum <- function(
+    x,
+    w,
+    hm,
+    n = TRUE,
+    h = TRUE,
+    env,
+    ...
+    ) {
+  
   # increment
   i <- get("i", envir = env)
 
@@ -77,10 +82,10 @@ rspectrum <- function(x,
     # in accordance to ploton et al. 2012
     output <- get("output", envir = env)
 
-    if (w / 2 < 29) {
-      output[i, ] <- rspec[1:floor(w / 2)]
+    if (length(rspec) < hm) {
+      output[i, 1:length(rspec)] <- rspec
     } else {
-      output[i, ] <- rspec[1:29]
+      output[i,] <- rspec[1:hm]
     }
 
     # assign values to matrix
diff --git a/man/foto.Rd b/man/foto.Rd
index d68360c..532c01a 100644
--- a/man/foto.Rd
+++ b/man/foto.Rd
@@ -7,6 +7,7 @@
 foto(
   x,
   window_size = 61,
+  harmonics = 29,
   method = "zones",
   norm_spec = FALSE,
   high_pass = TRUE,
@@ -20,6 +21,8 @@ foto(
 
 \item{window_size}{a moving window size in pixels (default = 61 pixels)}
 
+\item{harmonics}{number of harmonics to consider (29 by default)}
+
 \item{method}{zones (for discrete zones) or mw for a moving window
 approach}
 
diff --git a/man/rspectrum.Rd b/man/rspectrum.Rd
index cc864c1..019abe0 100644
--- a/man/rspectrum.Rd
+++ b/man/rspectrum.Rd
@@ -4,13 +4,15 @@
 \alias{rspectrum}
 \title{Calculates a radial spectrum}
 \usage{
-rspectrum(x, w, n = TRUE, h = TRUE, env, ...)
+rspectrum(x, w, hm, n = TRUE, h = TRUE, env, ...)
 }
 \arguments{
 \item{x}{a square matrix}
 
 \item{w}{a moving window size}
 
+\item{hm}{harmonics to consider as features}
+
 \item{n}{normalize, bolean \code{TRUE} or \code{FALSE}}
 
 \item{h}{high pass filter on the two first spectra values
diff --git a/vignettes/foto-vignette.Rmd b/vignettes/foto-vignette.Rmd
index 7ff9a47..f68805c 100644
--- a/vignettes/foto-vignette.Rmd
+++ b/vignettes/foto-vignette.Rmd
@@ -91,19 +91,18 @@ output_2 <- foto(r_2,
 ```
 
 
-
 ```{r figure_3, fig.width=7, fig.height=4}
 par(mfrow = c(1,2))
 plot(r_2, col = grey.colors(100), axes = FALSE, legend = FALSE, box = FALSE)
 plot(r, col = grey.colors(100), axes = FALSE, legend = FALSE, box = FALSE)
 
 par(mfrow = c(1,2))
-plotRGB(output_batch[[1]], stretch = "hist")
-plotRGB(output_batch[[2]], stretch = "hist")
+plotRGB(output_batch[[1]], stretch = "hist", smooth = FALSE)
+plotRGB(output_batch[[2]], stretch = "hist", smooth = FALSE)
 
 par(mfrow = c(1,2))
-plotRGB(output_2$rgb, stretch = "hist")
-plotRGB(output$rgb, stretch = "hist")
+plotRGB(output_2$rgb, stretch = "hist", smooth = FALSE)
+plotRGB(output$rgb, stretch = "hist", smooth = FALSE)
 ```