diff --git a/DESCRIPTION b/DESCRIPTION
index fed36a23..1539140b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -65,6 +65,7 @@ Collate:
     'hxsurf.R'
     'im3d.R'
     'interactive.R'
+    'morphometry.R'
     'nat-data.R'
     'nat-package.R'
     'ndigest.R'
diff --git a/NAMESPACE b/NAMESPACE
index 4f521cc8..996b07da 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -132,6 +132,7 @@ S3method(prune_online,neuronlist)
 S3method(remotesync,default)
 S3method(remotesync,neuronlistfh)
 S3method(resample,neuron)
+S3method(resample,neuronlist)
 S3method(rootpoints,default)
 S3method(rootpoints,igraph)
 S3method(rootpoints,neuron)
@@ -139,6 +140,8 @@ S3method(scale,dotprops)
 S3method(scale,neuron)
 S3method(setdiff,default)
 S3method(setdiff,neuronlist)
+S3method(sholl_analysis,neuron)
+S3method(sholl_analysis,neuronlist)
 S3method(subset,dotprops)
 S3method(subset,hxsurf)
 S3method(subset,neuron)
@@ -268,6 +271,8 @@ export(nrrd.voxdims)
 export(nvertices)
 export(nview3d)
 export(origin)
+export(overlap)
+export(overlap_score)
 export(plane_coefficients)
 export(pointsinside)
 export(potential_synapses)
@@ -308,6 +313,7 @@ export(seglist2swc)
 export(segmentgraph)
 export(select_points)
 export(setdiff)
+export(sholl_analysis)
 export(simplify_reglist)
 export(smooth_neuron)
 export(spine)
diff --git a/R/morphometry.R b/R/morphometry.R
new file mode 100644
index 00000000..baf683d2
--- /dev/null
+++ b/R/morphometry.R
@@ -0,0 +1,99 @@
+# Functions for morphological analysis
+
+
+#' Generate a connectivity matrix based on euclidean distance between points
+#'
+#' @description Generates an 'overlap matrix' of overlap scores between neurons in the \code{outputneurons} and \code{inputneurons} pools.
+#' For every point in a given neuron in \code{outputneurons}, a distance score is calculated to every point in a neuron in \code{inputneurons}.
+#' The sum of this score is added to the final output matrix. The score is calculated as \code{e(-d^2/(2*delta^2))}, where d is the euclidean distance between the two points,
+#' and delta is the expected distance in um that is considered 'close'. It is recommended that the user resamples neurons before use, using \code{\link{resample}}.
+#'
+#' @param outputneurons first set of neurons
+#' @param inputneurons second set of neurons
+#' @param delta the distance (in um) at which a synapse might occur
+#' @param progress whether or not to have a progress bar
+#' 
+#' @examples
+#' \dontrun{ 
+#' # Calculate how much some neurons overlap with one another
+#' ## Example requires the package nat.flybrains
+#' Cell07PNs_overlap = overlap_score(outputneurons = Cell07PNs, inputneurons = Cell07PNs)
+#' 
+#' ## Plot the results
+#' heatmap(Cell07PNs_overlap)
+#' } 
+#' @return a matrix of overlap scores
+#' @seealso \code{\link{potential_synapses}}, \code{\link{resample}}
+#' @export
+overlap_score <- function(outputneurons, inputneurons, 
+                    delta =1, progress = TRUE){
+  score.matrix = matrix(0,nrow = length(outputneurons),ncol = length(inputneurons))
+  rownames(score.matrix) = names(outputneurons)
+  colnames(score.matrix) = names(inputneurons)
+  for (n in 1:length(outputneurons)){
+    a = xyzmatrix(outputneurons[[n]])
+    inputneurons.d = nlapply(inputneurons, xyzmatrix, .progress = "none")
+    s = sapply(inputneurons.d, function(x)sum(exp(-nabor::knn(query = a, data = x,k=nrow(x))$nn.dists^2/(2*delta^2)))) # Score similar to that in Schlegel et al. 2015
+    score.matrix[n,] = s
+    if(progress){
+      nat_progress(x = n/length(outputneurons)*100, message = "calculating overlap")
+    } 
+  }
+  score.matrix
+}
+
+
+#' Perform a sholl analysis on neuron skeletons
+#'
+#' @description Functions for Sholl analysis of neuronal skeletons
+#'
+#' @param x a neuron or neuronlist object
+#' @param start the origin from which spheres are grown for the Sholl analysis
+#' @param starting.radius the radius of the first sphere. Defaults to the radius step
+#' @param ending.radius the radius of the last sphere. If NULL the distance to the furthest dendritic point from the start point is taken
+#' @param radius.step the change in radius between successive spheres. Defaults to one 100th of the radius of the ending sphere
+#' @return a data.frame of spheres radii and the number of dendritic intersections at each radius
+#' @examples
+#' \dontrun{ 
+#' # Calculate how much some neurons overlap with one another
+#' ## Example requires the package nat.flybrains
+#' Cell07PNs_sholl = sholl_analysis(x = Cell07PNs, radius.step = 1, ending.radius = 100)
+#' head(Cell07PNs_sholl[[1]])
+#' } 
+#' @export
+#' @rdname sholl_analysis
+sholl_analysis <- function(x, start = colMeans(xyzmatrix(x)), 
+                           starting.radius = radius.step, ending.radius = 1000, 
+                           radius.step = ending.radius/100) UseMethod("sholl_analysis")
+#' @export
+#' @rdname sholl_analysis
+sholl_analysis.neuron <- function(x, start = colMeans(xyzmatrix(x)), 
+                           starting.radius = radius.step, ending.radius = 1000, 
+                           radius.step = ending.radius/100){
+  unit.vector <- function(x) {x / sqrt(sum(x^2))}
+  dend = x$d
+  dend$dists = nabor::knn(data = matrix(start,ncol=3), query = nat::xyzmatrix(x),k=1)$nn.dists
+  if(is.null(ending.radius)){
+    ending.radius = max(dend$dists)
+  }
+  radii = seq(from = starting.radius, to = ending.radius, by = radius.step)
+  sholl = data.frame(radii = radii, intersections = 0)
+  for(n in 1:length(radii)){
+    r = radii[n]
+    segments = x$SegList
+    for(segment in segments){
+      p = dend[segment,]
+      dists = (nabor::knn(data = matrix(start,ncol=3), query = nat::xyzmatrix(p),k=1)$nn.dists - r) >= 0
+      sholl[n,]$intersections = sholl[n,]$intersections + lengths(regmatches(paste(dists,collapse=""), gregexpr("TRUEFALSE|FALSETRUE", paste(dists,collapse=""))))
+    }
+  }
+  sholl
+}
+#' @export
+#' @rdname sholl_analysis
+sholl_analysis.neuronlist <- function(x, start = colMeans(xyzmatrix(x)), 
+                                  starting.radius = radius.step, ending.radius = 1000, 
+                                  radius.step = ending.radius/100){
+  nlapply(x, sholl_analysis.neuron, 
+          start = start, starting.radius = starting.radius, ending.radius = ending.radius, radius.step = radius.step)
+}
diff --git a/R/neuron.R b/R/neuron.R
index a9c32545..7d64f676 100644
--- a/R/neuron.R
+++ b/R/neuron.R
@@ -582,6 +582,12 @@ resample.neuron<-function(x, stepsize, ...) {
   as.neuron(swc, origin=match(x$StartPoint, old_ids))
 }
 
+#' @export
+#' @rdname resample
+resample.neuronlist<-function(x, stepsize, ...){
+  nlapply(x, resample, stepsize=stepsize, ...)
+}
+
 # Interpolate ordered 3D points (optionally w diameter)
 # NB returns NULL if unchanged (when too short or <=2 points) 
 # and only returns _internal_ points, omitting the head and tail of a segment
diff --git a/R/ngraph.R b/R/ngraph.R
index b1fade48..8fd92ae5 100644
--- a/R/ngraph.R
+++ b/R/ngraph.R
@@ -620,7 +620,6 @@ EdgeListFromSegList<-function(SegList){
 #' LH_arbour = prune_in_volume(x = Cell07PNs, surf = nat.flybrains::IS2NP.surf,
 #'   neuropil = "LH_L", OmitFailures = TRUE)
 #' }
-#' @inheritParams prune
 #' @return A pruned neuron/neuronlist object
 #' @seealso \code{\link{as.neuron.ngraph}}, \code{\link{subset.neuron}},
 #'   \code{\link{prune.neuron}}, \code{\link{prune}}
diff --git a/R/potential_synapses.R b/R/potential_synapses.R
index fe0ea3da..f47eae58 100644
--- a/R/potential_synapses.R
+++ b/R/potential_synapses.R
@@ -291,3 +291,50 @@ restrictToBounds<-function(a, bounds){
         a[,2]>=bounds[3] & a[,2]<=bounds[4] &
         a[,3]>=bounds[5] & a[,3]<=bounds[6], ]
 }
+
+
+#' Generate a connectivity matrix based on euclidean distance between points
+#'
+#' @description Generates an 'overlap matrix' of overlap scores between neurons
+#'   in the \code{output.neurons} and \code{input.neurons} pools. For every
+#'   point in a given neuron in \code{output.neurons}, a distance score is
+#'   calculated to every point in a neuron in \code{input.neurons}. The sum of
+#'   this score is added to the final output matrix. The score is calculated as
+#'   \code{e(-d^2/(2*delta^2))}, where d is the euclidean distance between the
+#'   two points, and delta is the expected distance in um that is considered
+#'   'close'. It is recommended that the user resamples neurons before use,
+#'   using \code{\link{resample}}.
+#'
+#' @param output.neurons first set of neurons
+#' @param input.neurons second set of neurons
+#' @param delta the distance (in um) at which a synapse might occur
+#' @param progress whether or not to have a progress bar
+#'
+#' @examples
+#' \dontrun{
+#' # Calculate how much some neurons overlap with one another
+#' ## Example requires the package nat.flybrains
+#' Cell07PNs_overlap = overlap(output.neurons = Cell07PNs, input.neurons = Cell07PNs)
+#'
+#' ## Plot the results
+#' heatmap(Cell07PNs_overlap)
+#' }
+#' @return a matrix of overlap scores
+#' @seealso \code{\link{potential_synapses}}, \code{\link{resample}}
+#' @export
+overlap <- function(output.neurons, input.neurons, delta =1, progress = TRUE){
+  score.matrix = matrix(0,nrow = length(output.neurons),ncol = length(input.neurons))
+  rownames(score.matrix) = names(output.neurons)
+  colnames(score.matrix) = names(input.neurons)
+  for (n in 1:length(output.neurons)){
+    a = xyzmatrix(output.neurons[[n]])
+    input.neurons.d = nlapply(input.neurons, xyzmatrix, .progress = "none")
+    s = sapply(input.neurons.d, function(x)sum(exp(-nabor::knn(query = a, data = x,k=nrow(x))$nn.dists^2/(2*delta^2)))) # Score similar to that in Schlegel et al. 2015
+    score.matrix[n,] = s
+    if(progress){
+      nat_progress(x = n/length(output.neurons)*100, message = "calculating overlap")
+    } 
+  }
+  score.matrix
+}
+
diff --git a/R/utils.R b/R/utils.R
index e396415d..727bdcd8 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -6,3 +6,13 @@ file.exists <- function(...) {
   }
   base::file.exists(x)
 }
+
+# hidden
+nat_progress <- function (x, max = 100, message = NULL) {
+  percent <- x / max * 100
+  cat(sprintf('\r|%-50s| ~%d%% %s',
+              paste(rep('=', percent / 2), collapse = ''),
+              floor(percent), message))
+  if (x == max)
+    cat('\n')
+}
diff --git a/man/overlap.Rd b/man/overlap.Rd
new file mode 100644
index 00000000..72947db3
--- /dev/null
+++ b/man/overlap.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/potential_synapses.R
+\name{overlap}
+\alias{overlap}
+\title{Generate a connectivity matrix based on euclidean distance between points}
+\usage{
+overlap(output.neurons, input.neurons, delta = 1, progress = TRUE)
+}
+\arguments{
+\item{output.neurons}{first set of neurons}
+
+\item{input.neurons}{second set of neurons}
+
+\item{delta}{the distance (in um) at which a synapse might occur}
+
+\item{progress}{whether or not to have a progress bar}
+}
+\value{
+a matrix of overlap scores
+}
+\description{
+Generates an 'overlap matrix' of overlap scores between neurons
+  in the \code{output.neurons} and \code{input.neurons} pools. For every
+  point in a given neuron in \code{output.neurons}, a distance score is
+  calculated to every point in a neuron in \code{input.neurons}. The sum of
+  this score is added to the final output matrix. The score is calculated as
+  \code{e(-d^2/(2*delta^2))}, where d is the euclidean distance between the
+  two points, and delta is the expected distance in um that is considered
+  'close'. It is recommended that the user resamples neurons before use,
+  using \code{\link{resample}}.
+}
+\examples{
+\dontrun{
+# Calculate how much some neurons overlap with one another
+## Example requires the package nat.flybrains
+Cell07PNs_overlap = overlap(output.neurons = Cell07PNs, input.neurons = Cell07PNs)
+
+## Plot the results
+heatmap(Cell07PNs_overlap)
+}
+}
+\seealso{
+\code{\link{potential_synapses}}, \code{\link{resample}}
+}
diff --git a/man/overlap_score.Rd b/man/overlap_score.Rd
new file mode 100644
index 00000000..b0f5e519
--- /dev/null
+++ b/man/overlap_score.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/morphometry.R
+\name{overlap_score}
+\alias{overlap_score}
+\title{Generate a connectivity matrix based on euclidean distance between points}
+\usage{
+overlap_score(outputneurons, inputneurons, delta = 1, progress = TRUE)
+}
+\arguments{
+\item{outputneurons}{first set of neurons}
+
+\item{inputneurons}{second set of neurons}
+
+\item{delta}{the distance (in um) at which a synapse might occur}
+
+\item{progress}{whether or not to have a progress bar}
+}
+\value{
+a matrix of overlap scores
+}
+\description{
+Generates an 'overlap matrix' of overlap scores between neurons in the \code{outputneurons} and \code{inputneurons} pools.
+For every point in a given neuron in \code{outputneurons}, a distance score is calculated to every point in a neuron in \code{inputneurons}.
+The sum of this score is added to the final output matrix. The score is calculated as \code{e(-d^2/(2*delta^2))}, where d is the euclidean distance between the two points,
+and delta is the expected distance in um that is considered 'close'. It is recommended that the user resamples neurons before use, using \code{\link{resample}}.
+}
+\examples{
+\dontrun{ 
+# Calculate how much some neurons overlap with one another
+## Example requires the package nat.flybrains
+Cell07PNs_overlap = overlap_score(outputneurons = Cell07PNs, inputneurons = Cell07PNs)
+
+## Plot the results
+heatmap(Cell07PNs_overlap)
+} 
+}
+\seealso{
+\code{\link{potential_synapses}}, \code{\link{resample}}
+}
diff --git a/man/resample.Rd b/man/resample.Rd
index 385371ec..a72a66c9 100644
--- a/man/resample.Rd
+++ b/man/resample.Rd
@@ -3,11 +3,14 @@
 \name{resample}
 \alias{resample}
 \alias{resample.neuron}
+\alias{resample.neuronlist}
 \title{Resample an object with a new spacing}
 \usage{
 resample(x, ...)
 
 \method{resample}{neuron}(x, stepsize, ...)
+
+\method{resample}{neuronlist}(x, stepsize, ...)
 }
 \arguments{
 \item{x}{An object to resample}
diff --git a/man/sholl_analysis.Rd b/man/sholl_analysis.Rd
new file mode 100644
index 00000000..98aca6f3
--- /dev/null
+++ b/man/sholl_analysis.Rd
@@ -0,0 +1,57 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/morphometry.R
+\name{sholl_analysis}
+\alias{sholl_analysis}
+\alias{sholl_analysis.neuron}
+\alias{sholl_analysis.neuronlist}
+\title{Perform a sholl analysis on neuron skeletons}
+\usage{
+sholl_analysis(
+  x,
+  start = colMeans(xyzmatrix(x)),
+  starting.radius = radius.step,
+  ending.radius = 1000,
+  radius.step = ending.radius/100
+)
+
+\method{sholl_analysis}{neuron}(
+  x,
+  start = colMeans(xyzmatrix(x)),
+  starting.radius = radius.step,
+  ending.radius = 1000,
+  radius.step = ending.radius/100
+)
+
+\method{sholl_analysis}{neuronlist}(
+  x,
+  start = colMeans(xyzmatrix(x)),
+  starting.radius = radius.step,
+  ending.radius = 1000,
+  radius.step = ending.radius/100
+)
+}
+\arguments{
+\item{x}{a neuron or neuronlist object}
+
+\item{start}{the origin from which spheres are grown for the Sholl analysis}
+
+\item{starting.radius}{the radius of the first sphere. Defaults to the radius step}
+
+\item{ending.radius}{the radius of the last sphere. If NULL the distance to the furthest dendritic point from the start point is taken}
+
+\item{radius.step}{the change in radius between successive spheres. Defaults to one 100th of the radius of the ending sphere}
+}
+\value{
+a data.frame of spheres radii and the number of dendritic intersections at each radius
+}
+\description{
+Functions for Sholl analysis of neuronal skeletons
+}
+\examples{
+\dontrun{ 
+# Calculate how much some neurons overlap with one another
+## Example requires the package nat.flybrains
+Cell07PNs_sholl = sholl_analysis(x = Cell07PNs, radius.step = 1, ending.radius = 100)
+head(Cell07PNs_sholl[[1]])
+} 
+}