diff --git a/R/data_input.R b/R/data_input.R
index 96b9fa36ac59ef8791a385b6aaf320d30f8743cb..18ded0a6f303a1dd370063f8bd4396b003eae2b7 100644
--- a/R/data_input.R
+++ b/R/data_input.R
@@ -22,6 +22,7 @@ locs_from_csv <- function(file = NULL,
                           locprec.filter = 0,
                           locprecz.filter = 0) {
   # points <- read.csv(file, header = TRUE, stringsAsFactors = FALSE)
+  warning(paste0("Reading file ", file, "\n"))
   points <- data.table::fread(file, header = TRUE, stringsAsFactors = FALSE)
   colnames(points) <- sub("^xnm$", "x", colnames(points))
   colnames(points) <- sub("^ynm$", "y", colnames(points))
diff --git a/R/registration.R b/R/registration.R
index 75740f63f485ab4fad662895c6da496d6caddc4e..0fbbc0a5b344e5c89242e73113cf2c0c1917e604 100644
--- a/R/registration.R
+++ b/R/registration.R
@@ -810,9 +810,9 @@ jrmpc <-function(V, C = NULL, K = NULL, g = NULL, initialPriors = NULL, updatePr
   
   N <- sum(sapply(V, nrow))
   
-  ## Parameter h in the paper should be proportional to the volume that
+  ## Parameter h should be proportional to the volume that
   ## encompasses all the point sets.
-  ## Set to the volume of the enclosing sphere.
+  ## Set to it the volume of the enclosing sphere.
   if(d == 2) {
     h = pi*radius^2
   } else {
diff --git a/R/utils.R b/R/utils.R
index fa500be0970e42608f4855e8f76f578cd2a7d04c..58eac2ec11bf5b0566af85077c00a7fb447ce32f 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -325,7 +325,7 @@ points2img <- function(points, voxel.size, method, channels = NULL, ncpu = 1) {
       }
     } else {
       bins <- aws::binning(x = points[,c('x','y','z')], y = NULL, nbins = image.size)
-      I <- bins$table.freq
+      I <- array(bins$table.freq, dim = dim(bins$table.freq))
     }
   } else {
     stop(paste0("ERROR: Unknown method: ", method))
@@ -367,16 +367,47 @@ points_from_roi <- function(points, roi) {
 #' @param eps DBSCAN parameter, size of the epsilon neighbourhood
 #' @param minPts DBSCAN parameter, number of minimum points in the eps region
 #' @param keep.locprec logical (default: TRUE), whether to preserve the localization precision columns
+#' @param keep.channel logical (default: TRUE), whether to preserve channel information column
+#' @param cluster.2d logical (default: FALSE), whether to cluster only using x,y (and ignore z)
 #' @return a list of matrices with columns x,y,z and eventually locprec[z] and names set to the cluster indices.
 #' @export
 
-locs2ps <- function(points, eps, minPts, keep.locprec = TRUE) {
-  if(keep.locprec) {
-    points <- points[, c("x", "y", "z", colnames(points)[grep("locprec", colnames(points))])]
+locs2ps <- function(points, eps, minPts, keep.locprec = TRUE, keep.channel = TRUE, cluster.2d = FALSE) {
+  
+  is.2d <- FALSE
+  if(!("z" %in% colnames(points))) {
+    is.2d <- TRUE
+  }
+  if(keep.locprec && keep.channel) {
+    if(is.2d) {
+      points <- points[, c("x", "y", colnames(points)[grep("locprec|channel", colnames(points))])]
+    } else {
+      points <- points[, c("x", "y", "z", colnames(points)[grep("locprec|channel", colnames(points))])]
+    }
+  } else if(keep.locprec && !keep.channel) {
+    if(is.2d) {
+      points <- points[, c("x", "y", colnames(points)[grep("locprec", colnames(points))])]
+    } else {
+      points <- points[, c("x", "y", "z", colnames(points)[grep("locprec", colnames(points))])]
+    }
+  } else if(!keep.locprec && keep.channel) {
+    if(is.2d) {
+      points <- points[, c("x", "y", colnames(points)[grep("channel", colnames(points))])]
+    } else {
+      points <- points[, c("x", "y", "z", colnames(points)[grep("channel", colnames(points))])]
+    }
+  } else {
+    if(is.2d) {
+      points <- points[, c("x", "y")]
+    } else {
+      points <- points[, c("x", "y", "z")]
+    }
+  }
+  if(cluster.2d || is.2d) {
+    dbr <- dbscan::dbscan(points[, c("x", "y")], eps = eps, minPts = minPts, borderPoints = FALSE)
   } else {
-    points <- points[, c("x", "y", "z")]
+    dbr <- dbscan::dbscan(points[, c("x", "y", "z")], eps = eps, minPts = minPts, borderPoints = FALSE)
   }
-  dbr <- dbscan::dbscan(points[, c("x", "y", "z")], eps = eps, minPts = minPts, borderPoints = FALSE)
   points$site <- dbr$cluster
   idx.to.remove <- which(points$site==0, arr.ind = TRUE)
   if(length(idx.to.remove)>0) {
diff --git a/man/locs2ps.Rd b/man/locs2ps.Rd
index e18cf3a17632f2df1340a465c1ac0922d69cf74f..687de1a390f172a1760373ce44905c3b4d5ea111 100644
--- a/man/locs2ps.Rd
+++ b/man/locs2ps.Rd
@@ -4,7 +4,14 @@
 \alias{locs2ps}
 \title{locs2ps}
 \usage{
-locs2ps(points, eps, minPts, keep.locprec = TRUE)
+locs2ps(
+  points,
+  eps,
+  minPts,
+  keep.locprec = TRUE,
+  keep.channel = TRUE,
+  cluster.2d = FALSE
+)
 }
 \arguments{
 \item{points}{a point set as a data frame of coordinates with columns x,y,z.}
@@ -14,6 +21,10 @@ locs2ps(points, eps, minPts, keep.locprec = TRUE)
 \item{minPts}{DBSCAN parameter, number of minimum points in the eps region}
 
 \item{keep.locprec}{logical (default: TRUE), whether to preserve the localization precision columns}
+
+\item{keep.channel}{logical (default: TRUE), whether to preserve channel information column}
+
+\item{cluster.2d}{logical (default: FALSE), whether to cluster only using x,y (and ignore z)}
 }
 \value{
 a list of matrices with columns x,y,z and eventually locprec[z] and names set to the cluster indices.