diff --git a/.gitignore b/.gitignore
index e16642c..482c16f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -117,3 +117,7 @@ venv.bak/
 .idea/
 .vs/
 *.~lock.*
+
+
+# macOS
+*.DS_Store
diff --git a/doc_build/main/dev/CONTRIBUTING.md b/doc_build/main/dev/CONTRIBUTING.md
index eb8c61e..33795ff 100644
--- a/doc_build/main/dev/CONTRIBUTING.md
+++ b/doc_build/main/dev/CONTRIBUTING.md
@@ -28,7 +28,7 @@ clean commit history.
 
 
 #### Formatting
-`neuro` uses [Black](https://github.com/python/black) o ensure a consistent 
+`neuro` uses [Black](https://github.com/python/black) to ensure a consistent 
 code style. Please run `black ./ -l 79 --target-version py37` before making 
 any commits. To prevent any errors, it is easier to add a formatting check 
 as a [pre-commit hook](https://www.atlassian.com/git/tutorials/git-hooks). 
diff --git a/doc_build/main/user_guide/tools/cells_to_brainrender.md b/doc_build/main/user_guide/tools/cells_to_brainrender.md
index 4a976b3..8e23dac 100644
--- a/doc_build/main/user_guide/tools/cells_to_brainrender.md
+++ b/doc_build/main/user_guide/tools/cells_to_brainrender.md
@@ -5,7 +5,7 @@ To convert cell positions (e.g. in `cells_in_standard_space.xml`) to
 a format that can be used in 
 [BrainRender](https://github.com/BrancoLab/BrainRender).
 
-
+### Usage
 ```bash
     points_to_brainrender cells_in_standard_space.xml exported_cells.h5
 ```
@@ -16,7 +16,8 @@ Run `points_to_brainrender -h` to see all options.
 #### Positional
 * [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder) 
 cells file to be converted
-* Output filename. Should end with '.h5'
+* Output filename. Should end with '.h5'. If the containing directory doesn't 
+exist, it will be created.
 
 
 #### The following options may also need to be used:
diff --git a/doc_build/main/user_guide/tools/heatmap.md b/doc_build/main/user_guide/tools/heatmap.md
new file mode 100644
index 0000000..b36b990
--- /dev/null
+++ b/doc_build/main/user_guide/tools/heatmap.md
@@ -0,0 +1,47 @@
+# Heatmap generation
+
+To generate a heatmap of detected cells, which shows cell distributions in a 
+more intuitive way that showing individual cell positions:
+
+
+<img src="https://raw.githubusercontent.com/SainsburyWellcomeCentre/cellfinder/master/resources/heatmap.png" alt="heatmap" width="500"/>
+
+*Overlay on raw data and segentation from 
+[amap](https://github.com/SainsburyWellcomeCentre) added separately*
+
+
+### Usage
+```bash
+    heatmap cell_classification.xml heatmap.nii raw_data registered_atlas.nii -x 2 -y 2 -z 5
+```
+
+### Arguments
+Run `heatmap -h` to see all options.
+
+#### Positional arguments
+* [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder) 
+classified cells file (usually `cell_classification.xml`)
+* Output filename. Should end with '.nii'. If the containing directory doesn't 
+exist, it will be created.
+* Path to raw data (just a single channel). Used to find the shape of the 
+raw image that the detected cell positions are defined in.
+* Registered atlas file from [amap](https://github.com/SainsburyWellcomeCentre)
+ (typically run automatically in 
+ [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder)). File 
+ is usually `registered_atlas.nii`.
+ 
+
+#### Keyword arguments
+* `-x` or `--x-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the first dimension, specified in um.
+* `-y` or `--y-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the second dimension, specified in um.
+* `-z` or `--z-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the third dimension, specified in um. 
+
+
+#### The following options may also need to be used:
+* `--heatmap-bin` Heatmap bin size (um of each edge of histogram cube)
+* `--heatmap-smoothing` Gaussian smoothing sigma, in um.
+* `--no-mask-figs` Don't mask the figures (removing any areas outside the 
+brain, from e.g. smoothing)
diff --git a/docs/_sources/main/dev/CONTRIBUTING.md.txt b/docs/_sources/main/dev/CONTRIBUTING.md.txt
index eb8c61e..33795ff 100644
--- a/docs/_sources/main/dev/CONTRIBUTING.md.txt
+++ b/docs/_sources/main/dev/CONTRIBUTING.md.txt
@@ -28,7 +28,7 @@ clean commit history.
 
 
 #### Formatting
-`neuro` uses [Black](https://github.com/python/black) o ensure a consistent 
+`neuro` uses [Black](https://github.com/python/black) to ensure a consistent 
 code style. Please run `black ./ -l 79 --target-version py37` before making 
 any commits. To prevent any errors, it is easier to add a formatting check 
 as a [pre-commit hook](https://www.atlassian.com/git/tutorials/git-hooks). 
diff --git a/docs/_sources/main/user_guide/tools/cells_to_brainrender.md.txt b/docs/_sources/main/user_guide/tools/cells_to_brainrender.md.txt
index 4a976b3..8e23dac 100644
--- a/docs/_sources/main/user_guide/tools/cells_to_brainrender.md.txt
+++ b/docs/_sources/main/user_guide/tools/cells_to_brainrender.md.txt
@@ -5,7 +5,7 @@ To convert cell positions (e.g. in `cells_in_standard_space.xml`) to
 a format that can be used in 
 [BrainRender](https://github.com/BrancoLab/BrainRender).
 
-
+### Usage
 ```bash
     points_to_brainrender cells_in_standard_space.xml exported_cells.h5
 ```
@@ -16,7 +16,8 @@ Run `points_to_brainrender -h` to see all options.
 #### Positional
 * [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder) 
 cells file to be converted
-* Output filename. Should end with '.h5'
+* Output filename. Should end with '.h5'. If the containing directory doesn't 
+exist, it will be created.
 
 
 #### The following options may also need to be used:
diff --git a/docs/_sources/main/user_guide/tools/heatmap.md.txt b/docs/_sources/main/user_guide/tools/heatmap.md.txt
new file mode 100644
index 0000000..b36b990
--- /dev/null
+++ b/docs/_sources/main/user_guide/tools/heatmap.md.txt
@@ -0,0 +1,47 @@
+# Heatmap generation
+
+To generate a heatmap of detected cells, which shows cell distributions in a 
+more intuitive way that showing individual cell positions:
+
+
+<img src="https://raw.githubusercontent.com/SainsburyWellcomeCentre/cellfinder/master/resources/heatmap.png" alt="heatmap" width="500"/>
+
+*Overlay on raw data and segentation from 
+[amap](https://github.com/SainsburyWellcomeCentre) added separately*
+
+
+### Usage
+```bash
+    heatmap cell_classification.xml heatmap.nii raw_data registered_atlas.nii -x 2 -y 2 -z 5
+```
+
+### Arguments
+Run `heatmap -h` to see all options.
+
+#### Positional arguments
+* [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder) 
+classified cells file (usually `cell_classification.xml`)
+* Output filename. Should end with '.nii'. If the containing directory doesn't 
+exist, it will be created.
+* Path to raw data (just a single channel). Used to find the shape of the 
+raw image that the detected cell positions are defined in.
+* Registered atlas file from [amap](https://github.com/SainsburyWellcomeCentre)
+ (typically run automatically in 
+ [Cellfinder](https://github.com/SainsburyWellcomeCentre/cellfinder)). File 
+ is usually `registered_atlas.nii`.
+ 
+
+#### Keyword arguments
+* `-x` or `--x-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the first dimension, specified in um.
+* `-y` or `--y-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the second dimension, specified in um.
+* `-z` or `--z-pixel-size` Pixel spacing of the data that the cells are 
+defined in, in the third dimension, specified in um. 
+
+
+#### The following options may also need to be used:
+* `--heatmap-bin` Heatmap bin size (um of each edge of histogram cube)
+* `--heatmap-smoothing` Gaussian smoothing sigma, in um.
+* `--no-mask-figs` Don't mask the figures (removing any areas outside the 
+brain, from e.g. smoothing)
diff --git a/docs/genindex.html b/docs/genindex.html
index 0c05df6..5a93507 100644
--- a/docs/genindex.html
+++ b/docs/genindex.html
@@ -81,6 +81,7 @@
             
               <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/heatmap.html">Heatmap generation</a></li>
 </ul>
 <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/dev/CONTRIBUTING.html">Contributing</a></li>
diff --git a/docs/index.html b/docs/index.html
index 3a81562..802168a 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -81,6 +81,7 @@
             
               <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/heatmap.html">Heatmap generation</a></li>
 </ul>
 <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/dev/CONTRIBUTING.html">Contributing</a></li>
@@ -165,6 +166,7 @@ <h2>Tools<a class="headerlink" href="#tools" title="Permalink to this headline">
 <div class="toctree-wrapper compound">
 <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/heatmap.html">Heatmap generation</a></li>
 </ul>
 </div>
 </div>
diff --git a/docs/main/dev/CONTRIBUTING.html b/docs/main/dev/CONTRIBUTING.html
index 24e23a7..da05f65 100644
--- a/docs/main/dev/CONTRIBUTING.html
+++ b/docs/main/dev/CONTRIBUTING.html
@@ -35,7 +35,7 @@
   <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
     <link rel="index" title="Index" href="../../genindex.html" />
     <link rel="search" title="Search" href="../../search.html" />
-    <link rel="prev" title="Converting cells from cellfinder to brainrender format" href="../user_guide/tools/cells_to_brainrender.html" /> 
+    <link rel="prev" title="Heatmap generation" href="../user_guide/tools/heatmap.html" /> 
 </head>
 
 <body class="wy-body-for-nav">
@@ -81,6 +81,7 @@
             
               <ul>
 <li class="toctree-l1"><a class="reference internal" href="../user_guide/tools/cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1"><a class="reference internal" href="../user_guide/tools/heatmap.html">Heatmap generation</a></li>
 </ul>
 <ul class="current">
 <li class="toctree-l1 current"><a class="current reference internal" href="#">Contributing</a><ul>
@@ -186,7 +187,7 @@ <h2>Pull requests<a class="headerlink" href="#pull-requests" title="Permalink to
 </div>
 <div class="section" id="formatting">
 <h2>Formatting<a class="headerlink" href="#formatting" title="Permalink to this headline">¶</a></h2>
-<p><code class="docutils literal notranslate"><span class="pre">neuro</span></code> uses <a class="reference external" href="https://github.com/python/black">Black</a> o ensure a consistent
+<p><code class="docutils literal notranslate"><span class="pre">neuro</span></code> uses <a class="reference external" href="https://github.com/python/black">Black</a> to ensure a consistent
 code style. Please run <code class="docutils literal notranslate"><span class="pre">black</span> <span class="pre">./</span> <span class="pre">-l</span> <span class="pre">79</span> <span class="pre">--target-version</span> <span class="pre">py37</span></code> before making
 any commits. To prevent any errors, it is easier to add a formatting check
 as a <a class="reference external" href="https://www.atlassian.com/git/tutorials/git-hooks">pre-commit hook</a>.
@@ -235,7 +236,7 @@ <h2>Dependencies<a class="headerlink" href="#dependencies" title="Permalink to t
     <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
       
       
-        <a href="../user_guide/tools/cells_to_brainrender.html" class="btn btn-neutral float-left" title="Converting cells from cellfinder to brainrender format" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+        <a href="../user_guide/tools/heatmap.html" class="btn btn-neutral float-left" title="Heatmap generation" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
       
     </div>
   
diff --git a/docs/main/user_guide/tools/cells_to_brainrender.html b/docs/main/user_guide/tools/cells_to_brainrender.html
index de69fb3..bf15ea4 100644
--- a/docs/main/user_guide/tools/cells_to_brainrender.html
+++ b/docs/main/user_guide/tools/cells_to_brainrender.html
@@ -35,7 +35,7 @@
   <link rel="stylesheet" href="../../../_static/pygments.css" type="text/css" />
     <link rel="index" title="Index" href="../../../genindex.html" />
     <link rel="search" title="Search" href="../../../search.html" />
-    <link rel="next" title="Contributing" href="../../dev/CONTRIBUTING.html" />
+    <link rel="next" title="Heatmap generation" href="heatmap.html" />
     <link rel="prev" title="neuro" href="../../../index.html" /> 
 </head>
 
@@ -82,6 +82,7 @@
             
               <ul class="current">
 <li class="toctree-l1 current"><a class="current reference internal" href="#">Converting cells from cellfinder to brainrender format</a><ul>
+<li class="toctree-l2"><a class="reference internal" href="#usage">Usage</a></li>
 <li class="toctree-l2"><a class="reference internal" href="#arguments">Arguments</a><ul>
 <li class="toctree-l3"><a class="reference internal" href="#positional">Positional</a></li>
 <li class="toctree-l3"><a class="reference internal" href="#the-following-options-may-also-need-to-be-used">The following options may also need to be used:</a></li>
@@ -90,6 +91,7 @@
 <li class="toctree-l2"><a class="reference internal" href="#to-visualise-this-file-in-brainrender">To visualise this file in BrainRender</a></li>
 </ul>
 </li>
+<li class="toctree-l1"><a class="reference internal" href="heatmap.html">Heatmap generation</a></li>
 </ul>
 <ul>
 <li class="toctree-l1"><a class="reference internal" href="../../dev/CONTRIBUTING.html">Contributing</a></li>
@@ -160,9 +162,12 @@ <h1>Converting cells from cellfinder to brainrender format<a class="headerlink"
 <p>To convert cell positions (e.g. in <code class="docutils literal notranslate"><span class="pre">cells_in_standard_space.xml</span></code>) to
 a format that can be used in
 <a class="reference external" href="https://github.com/BrancoLab/BrainRender">BrainRender</a>.</p>
+<div class="section" id="usage">
+<h2>Usage<a class="headerlink" href="#usage" title="Permalink to this headline">¶</a></h2>
 <div class="highlight-bash notranslate"><div class="highlight"><pre><span></span>    points_to_brainrender cells_in_standard_space.xml exported_cells.h5
 </pre></div>
 </div>
+</div>
 <div class="section" id="arguments">
 <h2>Arguments<a class="headerlink" href="#arguments" title="Permalink to this headline">¶</a></h2>
 <p>Run <code class="docutils literal notranslate"><span class="pre">points_to_brainrender</span> <span class="pre">-h</span></code> to see all options.</p>
@@ -171,7 +176,8 @@ <h3>Positional<a class="headerlink" href="#positional" title="Permalink to this
 <ul class="simple">
 <li><p><a class="reference external" href="https://github.com/SainsburyWellcomeCentre/cellfinder">Cellfinder</a>
 cells file to be converted</p></li>
-<li><p>Output filename. Should end with ‘.h5’</p></li>
+<li><p>Output filename. Should end with ‘.h5’. If the containing directory doesn’t
+exist, it will be created.</p></li>
 </ul>
 </div>
 <div class="section" id="the-following-options-may-also-need-to-be-used">
@@ -208,7 +214,7 @@ <h2>To visualise this file in BrainRender<a class="headerlink" href="#to-visuali
   
     <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
       
-        <a href="../../dev/CONTRIBUTING.html" class="btn btn-neutral float-right" title="Contributing" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
+        <a href="heatmap.html" class="btn btn-neutral float-right" title="Heatmap generation" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
       
       
         <a href="../../../index.html" class="btn btn-neutral float-left" title="neuro" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
diff --git a/docs/main/user_guide/tools/heatmap.html b/docs/main/user_guide/tools/heatmap.html
new file mode 100644
index 0000000..56411ad
--- /dev/null
+++ b/docs/main/user_guide/tools/heatmap.html
@@ -0,0 +1,261 @@
+
+
+<!DOCTYPE html>
+<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
+<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
+<head>
+  <meta charset="utf-8">
+  
+  <meta name="viewport" content="width=device-width, initial-scale=1.0">
+  
+  <title>Heatmap generation &mdash; neuro  documentation</title>
+  
+
+  
+  
+  
+  
+
+  
+  <script type="text/javascript" src="../../../_static/js/modernizr.min.js"></script>
+  
+    
+      <script type="text/javascript" id="documentation_options" data-url_root="../../../" src="../../../_static/documentation_options.js"></script>
+        <script type="text/javascript" src="../../../_static/jquery.js"></script>
+        <script type="text/javascript" src="../../../_static/underscore.js"></script>
+        <script type="text/javascript" src="../../../_static/doctools.js"></script>
+        <script type="text/javascript" src="../../../_static/language_data.js"></script>
+    
+    <script type="text/javascript" src="../../../_static/js/theme.js"></script>
+
+    
+
+  
+  <link rel="stylesheet" href="../../../_static/css/theme.css" type="text/css" />
+  <link rel="stylesheet" href="../../../_static/pygments.css" type="text/css" />
+    <link rel="index" title="Index" href="../../../genindex.html" />
+    <link rel="search" title="Search" href="../../../search.html" />
+    <link rel="next" title="Contributing" href="../../dev/CONTRIBUTING.html" />
+    <link rel="prev" title="Converting cells from cellfinder to brainrender format" href="cells_to_brainrender.html" /> 
+</head>
+
+<body class="wy-body-for-nav">
+
+   
+  <div class="wy-grid-for-nav">
+    
+    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
+      <div class="wy-side-scroll">
+        <div class="wy-side-nav-search" >
+          
+
+          
+            <a href="../../../index.html" class="icon icon-home"> neuro
+          
+
+          
+          </a>
+
+          
+            
+            
+          
+
+          
+<div role="search">
+  <form id="rtd-search-form" class="wy-form" action="../../../search.html" method="get">
+    <input type="text" name="q" placeholder="Search docs" />
+    <input type="hidden" name="check_keywords" value="yes" />
+    <input type="hidden" name="area" value="default" />
+  </form>
+</div>
+
+          
+        </div>
+
+        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
+          
+            
+            
+              
+            
+            
+              <ul class="current">
+<li class="toctree-l1"><a class="reference internal" href="cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1 current"><a class="current reference internal" href="#">Heatmap generation</a><ul>
+<li class="toctree-l2"><a class="reference internal" href="#usage">Usage</a></li>
+<li class="toctree-l2"><a class="reference internal" href="#arguments">Arguments</a><ul>
+<li class="toctree-l3"><a class="reference internal" href="#positional-arguments">Positional arguments</a></li>
+<li class="toctree-l3"><a class="reference internal" href="#keyword-arguments">Keyword arguments</a></li>
+<li class="toctree-l3"><a class="reference internal" href="#the-following-options-may-also-need-to-be-used">The following options may also need to be used:</a></li>
+</ul>
+</li>
+</ul>
+</li>
+</ul>
+<ul>
+<li class="toctree-l1"><a class="reference internal" href="../../dev/CONTRIBUTING.html">Contributing</a></li>
+</ul>
+
+            
+          
+        </div>
+      </div>
+    </nav>
+
+    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">
+
+      
+      <nav class="wy-nav-top" aria-label="top navigation">
+        
+          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
+          <a href="../../../index.html">neuro</a>
+        
+      </nav>
+
+
+      <div class="wy-nav-content">
+        
+        <div class="rst-content">
+        
+          
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+<div role="navigation" aria-label="breadcrumbs navigation">
+
+  <ul class="wy-breadcrumbs">
+    
+      <li><a href="../../../index.html">Docs</a> &raquo;</li>
+        
+      <li>Heatmap generation</li>
+    
+    
+      <li class="wy-breadcrumbs-aside">
+        
+            
+        
+      </li>
+    
+  </ul>
+
+  
+  <hr/>
+</div>
+          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
+           <div itemprop="articleBody">
+            
+  <div class="section" id="heatmap-generation">
+<h1>Heatmap generation<a class="headerlink" href="#heatmap-generation" title="Permalink to this headline">¶</a></h1>
+<p>To generate a heatmap of detected cells, which shows cell distributions in a
+more intuitive way that showing individual cell positions:</p>
+<img src="https://raw.githubusercontent.com/SainsburyWellcomeCentre/cellfinder/master/resources/heatmap.png" alt="heatmap" width="500"/><p><em>Overlay on raw data and segentation from
+<a class="reference external" href="https://github.com/SainsburyWellcomeCentre">amap</a> added separately</em></p>
+<div class="section" id="usage">
+<h2>Usage<a class="headerlink" href="#usage" title="Permalink to this headline">¶</a></h2>
+<div class="highlight-bash notranslate"><div class="highlight"><pre><span></span>    heatmap cell_classification.xml heatmap.nii raw_data registered_atlas.nii -x <span class="m">2</span> -y <span class="m">2</span> -z <span class="m">5</span>
+</pre></div>
+</div>
+</div>
+<div class="section" id="arguments">
+<h2>Arguments<a class="headerlink" href="#arguments" title="Permalink to this headline">¶</a></h2>
+<p>Run <code class="docutils literal notranslate"><span class="pre">heatmap</span> <span class="pre">-h</span></code> to see all options.</p>
+<div class="section" id="positional-arguments">
+<h3>Positional arguments<a class="headerlink" href="#positional-arguments" title="Permalink to this headline">¶</a></h3>
+<ul class="simple">
+<li><p><a class="reference external" href="https://github.com/SainsburyWellcomeCentre/cellfinder">Cellfinder</a>
+classified cells file (usually <code class="docutils literal notranslate"><span class="pre">cell_classification.xml</span></code>)</p></li>
+<li><p>Output filename. Should end with ‘.nii’. If the containing directory doesn’t
+exist, it will be created.</p></li>
+<li><p>Path to raw data (just a single channel). Used to find the shape of the
+raw image that the detected cell positions are defined in.</p></li>
+<li><p>Registered atlas file from <a class="reference external" href="https://github.com/SainsburyWellcomeCentre">amap</a>
+(typically run automatically in
+<a class="reference external" href="https://github.com/SainsburyWellcomeCentre/cellfinder">Cellfinder</a>). File
+is usually <code class="docutils literal notranslate"><span class="pre">registered_atlas.nii</span></code>.</p></li>
+</ul>
+</div>
+<div class="section" id="keyword-arguments">
+<h3>Keyword arguments<a class="headerlink" href="#keyword-arguments" title="Permalink to this headline">¶</a></h3>
+<ul class="simple">
+<li><p><code class="docutils literal notranslate"><span class="pre">-x</span></code> or <code class="docutils literal notranslate"><span class="pre">--x-pixel-size</span></code> Pixel spacing of the data that the cells are
+defined in, in the first dimension, specified in um.</p></li>
+<li><p><code class="docutils literal notranslate"><span class="pre">-y</span></code> or <code class="docutils literal notranslate"><span class="pre">--y-pixel-size</span></code> Pixel spacing of the data that the cells are
+defined in, in the second dimension, specified in um.</p></li>
+<li><p><code class="docutils literal notranslate"><span class="pre">-z</span></code> or <code class="docutils literal notranslate"><span class="pre">--z-pixel-size</span></code> Pixel spacing of the data that the cells are
+defined in, in the third dimension, specified in um.</p></li>
+</ul>
+</div>
+<div class="section" id="the-following-options-may-also-need-to-be-used">
+<h3>The following options may also need to be used:<a class="headerlink" href="#the-following-options-may-also-need-to-be-used" title="Permalink to this headline">¶</a></h3>
+<ul class="simple">
+<li><p><code class="docutils literal notranslate"><span class="pre">--heatmap-bin</span></code> Heatmap bin size (um of each edge of histogram cube)</p></li>
+<li><p><code class="docutils literal notranslate"><span class="pre">--heatmap-smoothing</span></code> Gaussian smoothing sigma, in um.</p></li>
+<li><p><code class="docutils literal notranslate"><span class="pre">--no-mask-figs</span></code> Don’t mask the figures (removing any areas outside the
+brain, from e.g. smoothing)</p></li>
+</ul>
+</div>
+</div>
+</div>
+
+
+           </div>
+           
+          </div>
+          <footer>
+  
+    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
+      
+        <a href="../../dev/CONTRIBUTING.html" class="btn btn-neutral float-right" title="Contributing" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
+      
+      
+        <a href="cells_to_brainrender.html" class="btn btn-neutral float-left" title="Converting cells from cellfinder to brainrender format" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+      
+    </div>
+  
+
+  <hr/>
+
+  <div role="contentinfo">
+    <p>
+        &copy; Copyright 2020, Sainsbury Wellcome Centre
+
+    </p>
+  </div> 
+
+</footer>
+
+        </div>
+      </div>
+
+    </section>
+
+  </div>
+  
+
+
+  <script type="text/javascript">
+      jQuery(function () {
+          SphinxRtdTheme.Navigation.enable(true);
+      });
+  </script>
+
+  
+  
+    
+   
+
+</body>
+</html>
\ No newline at end of file
diff --git a/docs/objects.inv b/docs/objects.inv
index c6b92d8..7440bb3 100644
Binary files a/docs/objects.inv and b/docs/objects.inv differ
diff --git a/docs/search.html b/docs/search.html
index b0f4d8c..9782754 100644
--- a/docs/search.html
+++ b/docs/search.html
@@ -81,6 +81,7 @@
             
               <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/cells_to_brainrender.html">Converting cells from cellfinder to brainrender format</a></li>
+<li class="toctree-l1"><a class="reference internal" href="main/user_guide/tools/heatmap.html">Heatmap generation</a></li>
 </ul>
 <ul>
 <li class="toctree-l1"><a class="reference internal" href="main/dev/CONTRIBUTING.html">Contributing</a></li>
diff --git a/docs/searchindex.js b/docs/searchindex.js
index a52c11b..4a73c71 100644
--- a/docs/searchindex.js
+++ b/docs/searchindex.js
@@ -1 +1 @@
-Search.setIndex({docnames:["index","main/dev/CONTRIBUTING","main/user_guide/tools/cells_to_brainrender"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,"sphinx.ext.viewcode":1,sphinx:56},filenames:["index.rst","main/dev/CONTRIBUTING.md","main/user_guide/tools/cells_to_brainrender.md"],objects:{},objnames:{},objtypes:{},terms:{"20gb":[],"2gb":[],"case":1,"default":2,"function":1,"import":2,"new":1,"true":[1,2],"try":1,The:[],These:1,abl:[],activ:[],add:1,add_cells_from_fil:2,adding:1,advanc:[],all:[1,2],allen:[],allthough:[],also:[],altern:[],amap:[0,1],analysi:0,ani:1,anoth:1,appropri:[],approv:1,atla:2,atlas:[],automat:1,avail:[],avoid:[],bash:[],befor:1,being:[],below:[],black:1,brain:0,brainrend:0,branch:1,build:1,call:2,can:2,carri:[],cell:0,cellfind:[0,1],cellfinder_download:[],cells_in_standard_spac:2,centr:0,chang:[1,2],check:1,choos:[],clean:1,cli:1,clone:1,cluster:[],code:1,collect:0,com:1,command:1,commit:1,companion:1,compat:[],comput:[],conf:[],connect:[],consist:1,contribut:0,convert:0,correct:[],creat:[],cudatoolkit:[],current:[],data:[0,2],date:1,defin:2,design:1,dev:1,develop:1,dimens:2,disk:[],document:[],done:[],driver:[],easier:1,edit:1,either:1,elsewher:[],end:2,ensur:1,error:1,etc:[],everyth:[],exist:[],exit:1,exported_cel:2,extent:2,extract:1,fetch:1,file:[],filenam:2,first:2,fork:1,format:0,found:[],from:0,futur:[],git:1,github:1,has:2,have:[],hdf:2,henc:1,here:[],highli:[],histori:1,hook:1,how:[],http:1,identifi:2,imag:0,imlib:1,infer:[],initi:1,instal:1,instruct:[],internet:[],invers:1,isol:[],issu:1,jupyt:2,keep:1,kei:2,latest:[],learn:[],like:[],line:1,link:[],linux:1,machin:[],mai:1,main:1,maintain:1,make:1,master:1,max:2,maximum:2,merg:1,miniconda3:[],miniconda:[],minim:[],modifi:[],mous:[],must:2,name:[],necessari:1,need:[],neuro:1,onc:[],onli:[],origin:[],other:[],out:[],output:2,packag:1,particularli:0,path:[],pip:[0,1],pixel:2,pleas:1,points_to_brainrend:2,pre:1,prevent:1,previous:[],process:0,push:1,py37:1,pytest:1,python:[],recommend:[],refer:[],registr:[],relev:1,rememb:[],remot:1,render:2,repositori:1,run:[1,2],sainsburi:0,sainsburywellcomecentr:1,save:[],scene:2,second:2,see:2,set:1,should:[1,2],size:2,softwar:1,sourc:[],space:2,specifi:2,squash:1,still:[],style:1,submit:1,suppli:[],system:[],tag:1,target:1,tell:[],test:[],thi:1,third:2,time:[],tmp:[],tool:1,twice:[],ubuntu:[],unit_test:1,updat:[],upon:1,upstream:1,use:1,used:[],uses:1,using:[],version:1,via:1,want:[],wellcom:0,when:[],where:[],which:[],whole:0,wide:[],work:1,workflow:1,would:[],x86_64:[],xml:2,you:[],your:1,your_usernam:1},titles:["neuro","Contributing","Converting cells from cellfinder to brainrender format"],titleterms:{For:[],The:2,about:[],addit:[],also:2,argument:2,atla:[],brainrend:2,cell:2,cellfind:2,classif:[],conda:[],contribut:1,convert:2,cuda:[],cudnn:[],depend:1,develop:0,download:[],environ:[],file:2,follow:2,format:[1,2],from:2,gpu:[],guid:[],instal:0,introduct:0,mai:2,model:[],need:2,neuro:0,option:2,posit:2,pull:1,releas:1,request:1,requir:[],set:[],setup:1,support:[],test:1,thi:2,tool:0,train:[],travi:1,used:2,user:[],visualis:2}})
\ No newline at end of file
+Search.setIndex({docnames:["index","main/dev/CONTRIBUTING","main/user_guide/tools/cells_to_brainrender","main/user_guide/tools/heatmap"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,"sphinx.ext.viewcode":1,sphinx:56},filenames:["index.rst","main/dev/CONTRIBUTING.md","main/user_guide/tools/cells_to_brainrender.md","main/user_guide/tools/heatmap.md"],objects:{},objnames:{},objtypes:{},terms:{"case":1,"default":2,"function":1,"import":2,"new":1,"true":[1,2],"try":1,These:1,Used:3,add:1,add_cells_from_fil:2,added:3,adding:1,all:[1,2,3],amap:[0,1,3],analysi:0,ani:[1,3],anoth:1,approv:1,area:3,atla:[2,3],automat:[1,3],befor:1,bin:3,black:1,boundari:[],brain:[0,3],brainrend:0,branch:1,build:1,call:2,can:2,cell:[0,3],cell_classif:3,cellfind:[0,1,3],cells_in_standard_spac:2,centr:0,chang:[1,2],channel:3,check:1,classifi:3,clean:1,cli:1,clone:1,code:1,collect:0,com:1,command:1,commit:1,companion:1,consist:1,contain:[2,3],contribut:0,convert:0,creat:[2,3],cube:3,data:[0,2,3],date:1,defin:[2,3],design:1,detect:3,dev:1,develop:1,dimens:[2,3],directori:[2,3],distribut:3,doesn:[2,3],don:3,each:3,easier:1,edg:3,edit:1,either:1,end:[2,3],ensur:1,error:1,exist:[2,3],exit:1,exported_cel:2,extent:2,extract:1,fetch:1,fig:3,figur:3,file:3,filenam:[2,3],find:3,first:[2,3],fork:1,format:0,from:[0,3],gaussian:3,gener:0,git:1,github:1,has:2,hdf:2,heatmap:0,henc:1,histogram:3,histori:1,hook:1,http:1,identifi:2,imag:[0,3],imlib:1,individu:3,initi:1,instal:1,intuit:3,invers:1,issu:1,jupyt:2,just:3,keep:1,kei:2,line:1,linux:1,mai:1,main:1,maintain:1,make:1,mask:3,master:1,max:2,maximum:2,merg:1,more:3,must:2,necessari:1,neuro:1,nii:3,output:[2,3],outsid:3,overlai:3,packag:1,particularli:0,path:3,pip:[0,1],pixel:[2,3],pleas:1,points_to_brainrend:2,pre:1,prevent:1,process:0,push:1,py37:1,pytest:1,raw:3,raw_data:3,regist:3,registered_atla:3,relev:1,remot:1,remov:3,render:2,repositori:1,run:[1,2,3],sainsburi:0,sainsburywellcomecentr:1,scene:2,second:[2,3],see:[2,3],segent:3,separ:3,set:1,shape:3,should:[1,2,3],show:3,sigma:3,singl:3,size:[2,3],smooth:3,softwar:1,space:[2,3],specifi:[2,3],squash:1,style:1,submit:1,tag:1,target:1,thi:1,third:[2,3],tool:1,typic:3,unit_test:1,upon:1,upstream:1,use:1,uses:1,usual:3,version:1,via:1,wai:3,wellcom:0,which:3,whole:0,work:1,workflow:1,xml:[2,3],your:1,your_usernam:1},titles:["neuro","Contributing","Converting cells from cellfinder to brainrender format","Heatmap generation"],titleterms:{The:[2,3],also:[2,3],argument:[2,3],brainrend:2,cell:2,cellfind:2,contribut:1,convert:2,depend:1,develop:0,file:2,follow:[2,3],format:[1,2],from:2,gener:3,heatmap:3,instal:0,introduct:0,keyword:3,mai:[2,3],need:[2,3],neuro:0,option:[2,3],posit:[2,3],pull:1,releas:1,request:1,setup:1,test:1,thi:2,tool:0,travi:1,usag:[2,3],used:[2,3],visualis:2}})
\ No newline at end of file
diff --git a/neuro/atlas_tools/__init__.py b/neuro/atlas_tools/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/neuro/atlas_tools/custom_atlas_structures.py b/neuro/atlas_tools/custom_atlas_structures.py
new file mode 100644
index 0000000..a25be42
--- /dev/null
+++ b/neuro/atlas_tools/custom_atlas_structures.py
@@ -0,0 +1,163 @@
+import pandas as pd
+import numpy as np
+from brainio import brainio
+from skimage.filters import gaussian
+
+from neuro import brain_render_tools
+
+
+def load_atlas_structures_csv(path_to_structures_csv):
+    df = pd.read_csv(path_to_structures_csv)
+    return df
+
+
+def get_atlas_ids(df, parent_key):
+    ids = [int(id) for id in df[df["parent_id"] == parent_key]["id"]]
+    return ids
+
+
+def create_hierarchy_paths(df):
+    """
+    creates paths of id hierarchy to match custom atlas to allen atlas
+    :param df:
+    :return:
+    """
+    all_paths = {}
+    all_ids = df["id"]
+    all_ids = all_ids[all_ids.notnull()]
+    for id in all_ids:
+        path = get_structure_parents(df, id)
+        path = get_path_string_standard_fmt(path)
+        all_paths.setdefault(id, path)
+    return all_paths
+
+
+def get_structure_parents(df, k, parent=None, all_parents=None, root_id=997):
+    """
+    gets all parent structures of the given id
+
+    :param df:
+    :param k:
+    :param parent:
+    :param all_parents:
+    :param root_id:
+    :return:
+    """
+    if parent is None:
+        all_parents = []
+        parent = df[df["id"] == k]["parent_id"].values[0]
+        all_parents.append(parent)
+
+    if int(parent) != root_id:
+        parent = df[df["id"] == parent]["parent_id"].values[0]
+        all_parents.append(parent)
+        return get_structure_parents(df, k, parent, all_parents)
+    return all_parents[::-1]
+
+
+def get_path_string_standard_fmt(all_parent_ids):
+    all_parent_ids = [str(i) for i in all_parent_ids]
+    return "/".join(all_parent_ids)
+
+
+def add_to_df(df, df_dict):
+    df.insert(len(df.keys()), "structure_id_path", np.full(len(df), np.nan))
+    for k, v in df_dict.items():
+        loc = df[df["id"] == k].index[0]
+        df["structure_id_path"].iloc[loc] = v
+    return df
+
+
+def get_all_structure_children(df, k):
+    all_ids = df["id"]
+    all_ids = all_ids[all_ids.notnull()]
+    all_children = []
+    for id in all_ids:
+        all_parents = get_structure_parents(df, id)
+        if k in all_parents:
+            all_children.append(id)
+    if len(all_children) == 0:
+        all_children.append(k)
+    return all_children
+
+
+def render_all_subregions(
+    atlas_id,
+    out_dir,
+    atlas_path,
+    structures_csv_path,
+    smooth_threshold=0.4,
+    smooth_sigma=10,
+):
+    """
+    renders all children structures of a given id
+
+    :param atlas_id:
+    :param out_dir:
+    :param atlas_path:
+    :param structures_csv_path:
+    :return:
+    """
+    df = load_atlas_structures_csv(structures_csv_path)
+    atlas_ids = get_all_structure_children(df, atlas_id)
+    atlas = brainio.load_any(atlas_path)
+    for idx in atlas_ids:
+        region_mask = atlas == idx
+        smoothed_region = smooth_structure(
+            region_mask, threshold=smooth_threshold, sigma=smooth_sigma
+        )
+        brain_render_tools.volume_to_vector_array_to_obj_file(
+            smoothed_region, f"{out_dir}/{idx}.obj"
+        )
+
+
+def get_region_mask(atlas_id, atlas_path, structures_csv_path, smooth=False):
+    df = load_atlas_structures_csv(structures_csv_path)
+    atlas_ids = get_all_structure_children(df, atlas_id)
+    all_regions = get_region(atlas_ids, atlas_path, smooth=smooth)
+    return all_regions
+
+
+def get_region(atlas_ids, atlas_path, smooth=True):
+    atlas = brainio.load_any(atlas_path)
+    all_regions = np.zeros_like(atlas)
+    for id in atlas_ids:
+        region_mask = atlas == id
+        all_regions = np.logical_or(region_mask, all_regions)
+    if smooth:
+        all_regions = smooth_structure(all_regions)
+    return all_regions
+
+
+def get_arbitrary_structure_mask_from_custom_atlas(
+    atlas_ids, atlas_path, sigma, smoothing_threshold
+):
+    """
+    uses id numbers from custom atlas to generate a smoothed mask of those ids combined
+
+    :param atlas_ids:
+    :param atlas_path:
+    :param sigma:
+    :param smoothing_threshold:
+    :return:
+    """
+    atlas = brainio.load_any(atlas_path)
+    all_regions = np.zeros_like(atlas)
+    for id in atlas_ids:
+        region_mask = atlas == id
+        all_regions = np.logical_or(region_mask, all_regions)
+    all_regions = smooth_structure(
+        all_regions, threshold=smoothing_threshold, sigma=sigma
+    )
+    return all_regions
+
+
+def smooth_structure(image, threshold=0.4, sigma=10):
+    for i in range(image.shape[1]):
+        image[:, i, :] = gaussian(image[:, i, :], sigma) > threshold
+    return image
+
+
+def get_n_pixels_in_region(atlas_ids, atlas_path):
+    all_regions_smooth = get_region(atlas_ids, atlas_path)
+    return np.count_nonzero(all_regions_smooth)
diff --git a/neuro/atlas_tools/paths.py b/neuro/atlas_tools/paths.py
new file mode 100644
index 0000000..381ca4e
--- /dev/null
+++ b/neuro/atlas_tools/paths.py
@@ -0,0 +1,112 @@
+import os
+
+
+class Paths:
+    """
+    A single class to hold all file paths that amap may need. Any paths
+    prefixed with "tmp__" refer to internal intermediate steps, and will be
+    deleted if "--debug" is not used.
+    """
+
+    def __init__(self, registration_output_folder):
+        self.registration_output_folder = registration_output_folder
+        self.make_reg_paths()
+
+    def make_reg_paths(self):
+        self.brain_filtered = self.make_reg_path("brain_filtered.nii")
+
+        self.hemispheres = self.make_reg_path("hemispheres.nii")
+
+        self.annotations = self.make_reg_path("annotations.nii")
+
+        self.downsampled_brain_path = self.make_reg_path("downsampled.nii")
+        self.tmp__downsampled_filtered = self.make_reg_path(
+            "downsampled_filtered.nii"
+        )
+        self.registered_atlas_path = self.make_reg_path("registered_atlas.nii")
+        self.hemispheres_atlas_path = self.make_reg_path(
+            "registered_hemispheres.nii"
+        )
+        self.volume_csv_path = self.make_reg_path("volumes.csv")
+
+        self.tmp__affine_registered_atlas_brain_path = self.make_reg_path(
+            "affine_registered_atlas_brain.nii"
+        )
+        self.tmp__freeform_registered_atlas_brain_path = self.make_reg_path(
+            "freeform_registered_atlas_brain.nii"
+        )
+        self.tmp__inverse_freeform_registered_atlas_brain_path = self.make_reg_path(
+            "inverse_freeform_registered_brain.nii"
+        )
+
+        self.registered_atlas_img_path = self.make_reg_path(
+            "registered_atlas.nii"
+        )
+        self.registered_hemispheres_img_path = self.make_reg_path(
+            "registered_hemispheres.nii"
+        )
+
+        self.affine_matrix_path = self.make_reg_path("affine_matrix.txt")
+        self.invert_affine_matrix_path = self.make_reg_path(
+            "invert_affine_matrix.txt"
+        )
+
+        self.control_point_file_path = self.make_reg_path(
+            "control_point_file.nii"
+        )
+        self.inverse_control_point_file_path = self.make_reg_path(
+            "inverse_control_point_file.nii"
+        )
+
+        self.boundaries_file_path = self.make_reg_path("boundaries.nii")
+
+        (
+            self.tmp__affine_log_file_path,
+            self.tmp__affine_error_path,
+        ) = self.compute_reg_log_file_paths("affine")
+        (
+            self.tmp__freeform_log_file_path,
+            self.tmp__freeform_error_file_path,
+        ) = self.compute_reg_log_file_paths("freeform")
+        (
+            self.tmp__inverse_freeform_log_file_path,
+            self.tmp__inverse_freeform_error_file_path,
+        ) = self.compute_reg_log_file_paths("inverse_freeform")
+        (
+            self.tmp__segmentation_log_file,
+            self.tmp__segmentation_error_file,
+        ) = self.compute_reg_log_file_paths("segment")
+        (
+            self.tmp__invert_affine_log_file,
+            self.tmp__invert_affine_error_file,
+        ) = self.compute_reg_log_file_paths("invert_affine")
+
+    def make_reg_path(self, basename):
+        """
+        Compute the absolute path of the destination file to
+        self.registration_output_folder.
+
+        :param str basename:
+        :return: The path
+        :rtype: str
+        """
+        return os.path.join(self.registration_output_folder, basename)
+
+    def compute_reg_log_file_paths(self, basename):
+        """
+        Compute the path of the log and err file for the step corresponding
+        to basename
+
+        :param str basename:
+        :return: log_file_path, error_file_path
+        """
+
+        log_file_template = os.path.join(
+            self.registration_output_folder, "{}.log"
+        )
+        error_file_template = os.path.join(
+            self.registration_output_folder, "{}.err"
+        )
+        log_file_path = log_file_template.format(basename)
+        error_file_path = error_file_template.format(basename)
+        return log_file_path, error_file_path
diff --git a/neuro/brain_render_tools.py b/neuro/brain_render_tools.py
new file mode 100644
index 0000000..2b050c5
--- /dev/null
+++ b/neuro/brain_render_tools.py
@@ -0,0 +1,74 @@
+from brainrender.Utils.image import reorient_image, marching_cubes_to_obj
+from brainrender.scene import Scene
+from skimage import measure
+
+from neuro.atlas_tools.custom_atlas_structures import (
+    get_arbitrary_structure_mask_from_custom_atlas,
+)
+
+
+def render_region_from_custom_atlas(
+    output_dir,
+    atlas_ids,
+    structure_name,
+    atlas_path,
+    smoothing_threshold=0.4,
+    sigma=10,
+    voxel_size=10,
+):
+
+    all_regions = get_arbitrary_structure_mask_from_custom_atlas(
+        atlas_ids, atlas_path, sigma, smoothing_threshold
+    )
+
+    output_path = f"{output_dir}{structure_name}.obj"
+    oriented_binary = reorient_image(
+        all_regions, invert_axes=[2], orientation="coronal"
+    )
+
+    verts, faces, normals, values = measure.marching_cubes_lewiner(
+        oriented_binary, 0, step_size=1
+    )
+
+    if voxel_size is not 1:
+        verts = verts * voxel_size
+
+    faces = faces + 1
+    marching_cubes_to_obj((verts, faces, normals, values), output_path)
+
+
+def volume_to_vector_array_to_obj_file(image, output_path, voxel_size=10):
+
+    oriented_binary = reorient_image(
+        image, invert_axes=[2], orientation="coronal"
+    )
+
+    verts, faces, normals, values = measure.marching_cubes_lewiner(
+        oriented_binary, 0, step_size=1
+    )
+
+    if voxel_size is not 1:
+        verts = verts * voxel_size
+
+    faces = faces + 1
+    marching_cubes_to_obj((verts, faces, normals, values), str(output_path))
+
+
+def visualize_obj(obj_path, *args, color="lightcoral", **kwargs):
+    """
+        Uses brainrender to visualize a .obj file registered to the Allen CCF
+        :param obj_path: str, path to a .obj file
+        :param color: str, color of object being rendered
+    """
+    print("Visualizing : " + obj_path)
+    scene = Scene(add_root=True)
+    scene.add_from_file(obj_path, *args, c=color, **kwargs)
+
+    return scene
+
+
+def create_scene(default_structures):
+    scene = Scene(add_root=True)
+    for structure in default_structures:
+        scene.add_brain_regions([structure], use_original_color=True)
+    return scene
diff --git a/neuro/generic_neuro_tools.py b/neuro/generic_neuro_tools.py
new file mode 100644
index 0000000..101a650
--- /dev/null
+++ b/neuro/generic_neuro_tools.py
@@ -0,0 +1,139 @@
+import pathlib
+from pathlib import Path
+
+import numpy as np
+from brainio import brainio
+
+import imlib
+from cellfinder.summarise import count_summary as cells_regions
+from imlib.source.source_files import source_custom_config_cellfinder
+from imlib.general.exceptions import TransformationError
+from imlib.general.system import safe_execute_command, SafeExecuteCommandError
+from imlib.source.niftyreg_binaries import get_binary
+
+
+SOURCE_IMAGE_NAME = "downsampled.nii"
+DEFAULT_CONTROL_POINT_FILE = "inverse_control_point_file.nii"
+DEFAULT_OUTPUT_FILE_NAME = "roi_transformed.nii"
+DEFAULT_TEMP_FILE_NAME = "ROI_TMP.nii"
+PROGRAM_NAME = "reg_resample"
+
+
+def save_brain(image, source_image_path, output_path):
+    registration_config = source_custom_config_cellfinder()
+    atlas_scale, transformation_matrix = get_transform_space_params(
+        registration_config, source_image_path
+    )
+    brainio.to_nii(
+        image.astype(np.int16),
+        str(output_path),
+        scale=atlas_scale,
+        affine_transform=transformation_matrix,
+    )
+
+    print("finished generating roi image")
+
+
+def get_transform_space_params(registration_config, destination_image):
+    atlas = brainio.load_nii(str(destination_image), as_array=False)
+    atlas_scale = atlas.header.get_zooms()
+    atlas_pixel_sizes = cells_regions.get_atlas_pixel_sizes(
+        registration_config
+    )
+    transformation_matrix = np.eye(4)
+    for i, axis in enumerate(("x", "y", "z")):
+        transformation_matrix[i, i] = atlas_pixel_sizes[axis]
+    return atlas_scale, transformation_matrix
+
+
+def get_transformation_matrix(self):
+    atlas_pixel_sizes = cells_regions.get_atlas_pixel_sizes(self._atlas_config)
+    transformation_matrix = np.eye(4)
+    for i, axis in enumerate(("x", "y", "z")):
+        transformation_matrix[i, i] = atlas_pixel_sizes[axis]
+    self.transformation_matrix = transformation_matrix
+
+
+def get_registration_cmd(
+    program_path,
+    floating_image_path,
+    output_file_name,
+    destination_image_filename,
+    control_point_file,
+):
+    # TODO combine with amap.brain_registration
+    cmd = "{} -cpp {} -flo {} -ref {} -res {}".format(
+        program_path,
+        control_point_file,
+        floating_image_path,
+        destination_image_filename,
+        output_file_name,
+    )
+    return cmd
+
+
+def transform_image_to_standard_space(
+    reg_dir,
+    image_to_transform_fname="downsampled.nii",
+    output_fname="background_channel_reg_to_filtered_brain.nii",
+):
+
+    reg_dir = Path(reg_dir)
+    image_to_transform = reg_dir / image_to_transform_fname
+    image_with_desired_coordinate_space = reg_dir / "brain_filtered.nii"
+    control_point_file = reg_dir / "inverse_control_point_file.nii"
+    output_path = reg_dir / output_fname
+
+    nifty_reg_binaries_folder = (
+        imlib.source.niftyreg_binaries.get_niftyreg_binaries()
+    )
+    program_path = get_binary(nifty_reg_binaries_folder, PROGRAM_NAME)
+
+    reg_cmd = get_registration_cmd(
+        program_path,
+        image_to_transform,
+        output_path,
+        image_with_desired_coordinate_space,
+        control_point_file,
+    )
+
+    log_file_path = output_path.parent / "roi_transform_log.txt"
+    error_file_path = output_path.parent / "roi_transform_error.txt"
+    safely_execute_amap_registration(error_file_path, log_file_path, reg_cmd)
+    print(f"Registered ROI image can be found at {output_path}")
+
+
+def safely_execute_amap_registration(error_file_path, log_file_path, reg_cmd):
+    print("Running ROI registration")
+    try:
+        safe_execute_command(reg_cmd, log_file_path, error_file_path)
+    except SafeExecuteCommandError as err:
+        raise TransformationError("ROI registration failed; {}".format(err))
+
+
+def transform_all_channels_to_standard_space(reg_dir):
+    p = pathlib.Path(reg_dir)
+    all_channels = list(p.glob("downsampled_channel*"))
+    if any("registered" in f.name for f in all_channels):
+        return "looks like already processed... skipping..."
+
+    for channel in all_channels:
+        transform_image_to_standard_space(
+            reg_dir,
+            image_to_transform_fname=channel.name,
+            output_fname=f"registered_{channel.name}",
+        )
+
+
+def transform_background_channel_to_standard_space(reg_dir):
+    p = pathlib.Path(reg_dir)
+    all_channels = list(p.glob("downsampled.nii"))
+    if any("registered" in f.name for f in all_channels):
+        return "looks like already processed... skipping..."
+
+    for channel in all_channels:
+        transform_image_to_standard_space(
+            reg_dir,
+            image_to_transform_fname=channel.name,
+            output_fname=f"registered_{channel.name}",
+        )
diff --git a/neuro/heatmap/heatmap.py b/neuro/heatmap/heatmap.py
index 742ee62..2013df0 100644
--- a/neuro/heatmap/heatmap.py
+++ b/neuro/heatmap/heatmap.py
@@ -6,10 +6,9 @@
 import logging
 import argparse
 
+from pathlib import Path
 import numpy as np
 from skimage.filters import gaussian
-from skimage.transform import resize
-
 from brainio import brainio
 from imlib.cells.utils import get_cell_location_array
 from imlib.image.scale import scale_and_convert_to_16_bits
@@ -17,6 +16,8 @@
 from imlib.image.shape import convert_shape_dict_to_array_shape
 from imlib.image.masking import mask_image_threshold
 from imlib.general.numerical import check_positive_float
+from imlib.general.system import ensure_directory_exists
+from imlib.image.size import resize_array
 
 
 def run(
@@ -65,9 +66,12 @@ def run(
 
     logging.debug("Generating heatmap (3D histogram)")
     heatmap_array, _ = np.histogramdd(cells_array, bins=bins)
+    # otherwise resized array is too big to fit into RAM
+    heatmap_array = heatmap_array.astype(np.uint16)
 
     logging.debug("Resizing heatmap to the size of the target image")
-    heatmap_array = resize(heatmap_array, target_size, order=0)
+    heatmap_array = resize_array(heatmap_array, target_size)
+
     if smoothing is not None:
         logging.debug(
             "Applying Gaussian smoothing with a kernel sigma of: "
@@ -85,6 +89,9 @@ def run(
         logging.debug("Converting to 16 bit")
         heatmap_array = scale_and_convert_to_16_bits(heatmap_array)
 
+    logging.debug("Ensuring output directory exists")
+    ensure_directory_exists(Path(output_filename).parent)
+
     logging.debug("Saving heatmap image")
     brainio.to_nii(
         heatmap_array,
@@ -128,6 +135,7 @@ def get_parser():
         type=check_positive_float,
         help="Pixel spacing of the data in the first "
         "dimension, specified in um.",
+        required=True,
     )
     parser.add_argument(
         "-y",
@@ -136,6 +144,7 @@ def get_parser():
         type=check_positive_float,
         help="Pixel spacing of the data in the second "
         "dimension, specified in um.",
+        required=True,
     )
     parser.add_argument(
         "-z",
@@ -144,6 +153,7 @@ def get_parser():
         type=check_positive_float,
         help="Pixel spacing of the data in the third "
         "dimension, specified in um.",
+        required=True,
     )
     parser.add_argument(
         "--heatmap-smoothing",
@@ -230,7 +240,7 @@ def _get_binning(self):
         )
 
     def _get_atlas_data(self):
-        self.atlas_data = self._downsampled_image.get_fdata()
+        self.atlas_data = brainio.load_nii(self._target_image, as_array=True)
 
     def _get_atlas_scale(self):
         self.atlas_scale = self._downsampled_image.header.get_zooms()
diff --git a/neuro/injection_finder/__init__.py b/neuro/injection_finder/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/neuro/injection_finder/extraction.py b/neuro/injection_finder/extraction.py
new file mode 100644
index 0000000..29ac760
--- /dev/null
+++ b/neuro/injection_finder/extraction.py
@@ -0,0 +1,217 @@
+import os
+import numpy as np
+from pathlib import Path
+
+from skimage.filters import gaussian as gaussian_filter
+from skimage.filters import threshold_otsu
+from skimage import measure
+
+from brainio import brainio
+from imlib.IO.surfaces import marching_cubes_to_obj
+from imlib.image.orient import reorient_image
+from imlib.image.objects import keep_n_largest_objects
+
+from neuro.injection_finder.registration import get_registered_image
+from neuro.injection_finder.parsers import extraction_parser
+
+import neuro as package_for_log
+import logging
+from fancylog import fancylog
+
+
+class Extractor:
+    def __init__(
+        self,
+        img_filepath,
+        registration_folder,
+        overwrite=False,
+        gaussian_kernel=2,
+        gaussian_kernel_z=6,
+        percentile_threshold=99.95,
+        threshold_type="otsu",
+        obj_path=None,
+        overwrite_registration=False,
+    ):
+        """
+        Extractor processes a downsampled.nii image to extract the location of
+        the injection site.
+        This is done by registering the image to the allen CCF, blurring,
+        thresholding and finally a marching cube algorithm to extract the
+        surface of the injection site.
+
+        :param img_filepath: str, path to .nii file
+        :param registration_folder: str, path to the registration folder
+        [from cellfinder or amap]
+        :param overwrite: bool, if False it will avoid overwriting files
+        :param gaussian_kernel: float, size of kernel used for smoothing (x, y)
+        :param gaussian_kernel_z: float, size of kernel used for smoothing (z)
+        :param percentile_threshold: float, in range [0, 1] percentile to use
+        for thresholding
+        :param threshold_type: str, either ['otsu', 'percentile'],
+        type of threshold used
+        :param obj_path: path to .obj file destination.
+        :param overwrite_registration: if false doesn't overwrite the
+        registration step
+        """
+
+        # Get arguments
+        self.img_filepath = img_filepath
+        self.registration_folder = registration_folder
+        self.logging = logging
+        self.overwrite = overwrite
+        self.gaussian_kernel = gaussian_kernel
+        self.gaussian_kernel_z = gaussian_kernel_z
+        self.percentile_threshold = percentile_threshold
+        self.threshold_type = threshold_type
+        self.obj_path = obj_path
+        self.overwrite_registration = overwrite_registration
+
+        # Run first with the image oriented for brainrender
+        image = self.setup()
+        self.extract(image)
+
+    def setup(self):
+        if not os.path.isfile(self.img_filepath):
+            raise FileNotFoundError(
+                "The image path {} is not valid".format(self.img_filepath)
+            )
+
+        self.thresholded_savepath = (
+            self.img_filepath.split(".")[0] + "_thresholded.nii"
+        )
+
+        # Get path to obj file and check if it exists
+        if self.obj_path is None:
+            self.obj_path = self.img_filepath.split(".")[0] + ".obj"
+
+        if os.path.isfile(self.obj_path) and not self.overwrite:
+            self.logging.warning(
+                "A file exists already at {}."
+                "Analysis will not run as overwrite is set disabled".format(
+                    self.obj_path
+                )
+            )
+
+        # Load image and register
+        image = get_registered_image(
+            self.img_filepath,
+            self.registration_folder,
+            overwrite=self.overwrite_registration,
+        )
+        return image
+
+    def extract(self, image, voxel_size=10):
+        self.logging.info("Processing " + self.img_filepath)
+        self.logging.info(
+            "Gaussian filtering with kernel size: {}".format(
+                self.gaussian_kernel
+            )
+        )
+
+        # Gaussian filter
+        kernel_shape = [
+            self.gaussian_kernel,
+            self.gaussian_kernel,
+            self.gaussian_kernel_z,
+        ]
+        image = gaussian_filter(image, kernel_shape)
+        self.logging.info("Filtering completed")
+
+        # Thresholding
+        if self.threshold_type.lower() == "otsu":
+            thresh = threshold_otsu(image)
+            self.logging.info(
+                "Thresholding with {} threshold type".format(
+                    self.threshold_type
+                )
+            )
+
+        elif (
+            self.threshold_type.lower() == "percentile"
+            or self.threshold_type.lower() == "perc"
+        ):
+            thresh = np.percentile(image.ravel(), self.percentile_threshold)
+            self.logging.info(
+                "Thresholding with {} threshold type. "
+                "{}th percentile [{}]".format(
+                    self.threshold_type, self.percentile_threshold, thresh
+                )
+            )
+        else:
+            raise ValueError(
+                "Unrecognised thresholding type: " + self.threshold_type
+            )
+
+        binary = image > thresh
+        binary = keep_n_largest_objects(binary)
+
+        # Save thresholded image
+        if not os.path.isfile(self.thresholded_savepath) or self.overwrite:
+            self.logging.info(
+                "Saving thresholded image to {}".format(
+                    self.thresholded_savepath
+                )
+            )
+            brainio.to_nii(binary.astype(np.int16), self.thresholded_savepath)
+
+        binary = reorient_image(
+            binary, invert_axes=[2,], orientation="coronal"
+        )
+
+        # apply marching cubes
+        self.logging.info("Extracting surface from thresholded image")
+        verts, faces, normals, values = measure.marching_cubes_lewiner(
+            binary, 0, step_size=1
+        )
+
+        # Scale to atlas spacing
+        if voxel_size is not 1:
+            verts = verts * voxel_size
+
+        # Save image to .obj
+        self.logging.info(" Saving .obj at {}".format(self.obj_path))
+        faces = faces + 1
+        marching_cubes_to_obj((verts, faces, normals, values), self.obj_path)
+
+
+def main():
+    args = extraction_parser().parse_args()
+
+    # Get output directory
+    if args.output_directory is None:
+        outdir = os.getcwd()
+    elif not os.path.isdir(args.output_directory):
+        raise ValueError("Output directory invalid")
+    else:
+        outdir = args.output_directory
+
+    if args.obj_path is None:
+        args.obj_path = Path(args.img_filepath).with_suffix(".obj")
+    else:
+        args.obj_path = Path(args.obj_path)
+
+    # Start log
+    fancylog.start_logging(
+        outdir,
+        package_for_log,
+        filename="injection_finder",
+        verbose=args.debug,
+        log_to_file=args.save_log,
+    )
+
+    # Start extraction
+    Extractor(
+        args.img_filepath,
+        args.registration_folder,
+        overwrite=args.overwrite,
+        gaussian_kernel=args.gaussian_kernel,
+        gaussian_kernel_z=args.gaussian_kernel_z,
+        percentile_threshold=args.percentile_threshold,
+        threshold_type=args.threshold_type,
+        obj_path=args.obj_path,
+        overwrite_registration=args.overwrite_registration,
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/neuro/injection_finder/parsers.py b/neuro/injection_finder/parsers.py
new file mode 100644
index 0000000..2e3cba6
--- /dev/null
+++ b/neuro/injection_finder/parsers.py
@@ -0,0 +1,102 @@
+import argparse
+
+
+def extraction_parser():
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+    parser.add_argument(
+        dest="img_filepath", type=str, help="Path to brain volume (.nii) data",
+    )
+
+    parser.add_argument(
+        dest="registration_folder",
+        type=str,
+        help="Path to cellfinder registration folder",
+    )
+
+    parser.add_argument(
+        "-od",
+        "--output-directory",
+        dest="output_directory",
+        type=str,
+        default=None,
+        help="Path to directory where the log will be saved.",
+    )
+
+    parser.add_argument(
+        "-ow",
+        "--overwrite",
+        dest="overwrite",
+        type=bool,
+        default=True,
+        help="If false files will not be overwritten.",
+    )
+
+    parser.add_argument(
+        "-o",
+        "--obj-path",
+        dest="obj_path",
+        type=str,
+        default=None,
+        help="Path to output .obj file. Will default to the image directory.",
+    )
+
+    parser.add_argument(
+        "-k",
+        "--gaussian-kernel",
+        dest="gaussian_kernel",
+        type=float,
+        default=2,
+        help="Float, size of kernel for gaussian smoothing (x,y directions).",
+    )
+    parser.add_argument(
+        "--gaussian-kernel-z",
+        dest="gaussian_kernel_z",
+        type=float,
+        default=6,
+        help="Float, size of kernel for gaussian smoothing (zdirections).",
+    )
+
+    parser.add_argument(
+        "-pt",
+        "--percentile_threshold",
+        dest="percentile_threshold",
+        type=float,
+        default=99.995,
+        help="Float in range [0, 100]. The percentile number of pixel "
+        "intensity values for tresholding",
+    )
+
+    parser.add_argument(
+        "-tt",
+        "--treshold-type",
+        dest="threshold_type",
+        type=str,
+        default="otsu",
+        help="'otsu' or 'percentile'. Determines how the threshold "
+        "value is computed",
+    )
+
+    parser.add_argument(
+        "-or",
+        "--overwrite-registration",
+        dest="overwrite_registration",
+        type=str,
+        default="False",
+        help="If false skip running again the registration",
+    )
+    parser.add_argument(
+        "--debug",
+        dest="debug",
+        action="store_true",
+        help="Debug mode. Will increase verbosity of logging and save all "
+        "intermediate files for diagnosis of software issues.",
+    )
+    parser.add_argument(
+        "--save-log",
+        dest="save_log",
+        action="store_true",
+        help="Save logging to file (in addition to logging to terminal).",
+    )
+    return parser
diff --git a/neuro/injection_finder/registration.py b/neuro/injection_finder/registration.py
new file mode 100644
index 0000000..3abe057
--- /dev/null
+++ b/neuro/injection_finder/registration.py
@@ -0,0 +1,41 @@
+import os
+import logging
+
+from brainio import brainio
+from imlib.register.niftyreg.transform import run_transform
+
+
+DEFAULT_CONTROL_POINT_FILE = "inverse_control_point_file.nii"
+default_atlas_name = "brain_filtered.nii"
+
+
+def get_registered_image(nii_path, registration_dir, overwrite=False):
+    # get file paths
+    basedir = os.path.split(nii_path)[0]
+    output_filename = os.path.join(
+        basedir,
+        "{}_transformed.nii".format(os.path.split(nii_path)[1].split(".")[0]),
+    )
+    if os.path.isfile(output_filename) and not overwrite:
+        logging.info("Skipping registration as output file already exists")
+    else:
+
+        destination_image = os.path.join(registration_dir, default_atlas_name)
+        control_point_file = os.path.join(
+            registration_dir, DEFAULT_CONTROL_POINT_FILE
+        )
+
+        log_file_path = os.path.join(basedir, "registration_log.txt")
+        error_file_path = os.path.join(basedir, "registration_err.txt")
+
+        logging.info("Running registration")
+        run_transform(
+            nii_path,
+            output_filename,
+            destination_image,
+            control_point_file,
+            log_file_path,
+            error_file_path,
+        )
+
+    return brainio.load_any(output_filename)
diff --git a/neuro/lesion_and_track_tools/__init__.py b/neuro/lesion_and_track_tools/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/neuro/lesion_and_track_tools/fiber_tract_viewer.py b/neuro/lesion_and_track_tools/fiber_tract_viewer.py
new file mode 100644
index 0000000..943763d
--- /dev/null
+++ b/neuro/lesion_and_track_tools/fiber_tract_viewer.py
@@ -0,0 +1,90 @@
+import napari
+import pathlib
+import numpy as np
+
+from PySide2.QtWidgets import QApplication
+from neuro.visualise.vis_tools import prepare_load_nii
+from brainio import brainio
+
+from neuro.brain_render_tools import volume_to_vector_array_to_obj_file
+from neuro.lesion_and_track_tools.lesion_and_track_estimation import (
+    get_fiber_track,
+)
+from neuro.generic_neuro_tools import (
+    transform_background_channel_to_standard_space,
+)
+
+
+def display_channel(viewer, reg_dir, channel_fname):
+    """
+    Display downsampled data
+    :param viewer:
+    :param args:
+    :param paths:
+    :return:
+    """
+    reg_dir = pathlib.Path(reg_dir)
+
+    viewer.add_image(
+        prepare_load_nii(reg_dir / channel_fname, memory=False),
+        name="Downsampled filtered_brain_space",
+    )
+
+
+def run_track_viewer(
+    reg_dir,
+    atlas_name="annotations.nii",
+    background_channel_name="registered_downsampled.nii",
+    output_name="registered_track_{}.nii",
+):
+    reg_dir = pathlib.Path(reg_dir)
+    print("Starting amap viewer")
+
+    with napari.gui_qt():
+        v = napari.Viewer(title="track viewer")
+        display_channel(v, reg_dir, background_channel_name)
+        points_layer = v.add_points()
+
+        @v.bind_key("b")
+        def conv_to_b(v):
+            selected_points = v.layers[-1].selected_data
+            if len(selected_points) == 1:
+                print(f"labeled point at {v.layers[-1].data}")
+                z = int(v.layers[-1].data[0][0])
+                y = int(v.layers[-1].data[0][1])
+                x = int(v.layers[-1].data[0][2])
+                seed_point = (x, y, z)
+                seed_point_str = f"{x}_{y}_{z}"
+                print(f"{seed_point} extracted, getting fiber track from this")
+                get_fiber_track(
+                    reg_dir / atlas_name,
+                    reg_dir / background_channel_name,
+                    seed_point,
+                    reg_dir / output_name.format(seed_point_str),
+                )
+                print("closing application...")
+                QApplication.closeAllWindows()
+
+        @v.bind_key("q")
+        def quit(v):
+            QApplication.closeAllWindows()
+
+
+def load_arrays(images, names):
+    with napari.gui_qt():
+        v = napari.Viewer(title="track viewer")
+        for image, name in zip(images, names):
+            image = np.swapaxes(image, 2, 0)
+            v.add_image(image, name=name)
+
+
+def get_fiber_tract_in_standard_space(
+    reg_dir, output_name="registered_track.nii"
+):
+    transform_background_channel_to_standard_space(reg_dir)
+    run_track_viewer(reg_dir)
+    reg_dir = pathlib.Path(reg_dir)
+    brain = brainio.load_any(str(reg_dir / output_name))
+    volume_to_vector_array_to_obj_file(
+        brain, str(reg_dir / output_name.replace(".nii", ".obj"))
+    )
diff --git a/neuro/lesion_and_track_tools/lesion_and_track_estimation.py b/neuro/lesion_and_track_tools/lesion_and_track_estimation.py
new file mode 100644
index 0000000..c21427a
--- /dev/null
+++ b/neuro/lesion_and_track_tools/lesion_and_track_estimation.py
@@ -0,0 +1,118 @@
+import os
+import pathlib
+
+import numpy as np
+import skimage
+from skimage.segmentation import flood
+from skimage.morphology import binary_erosion, binary_dilation
+from skimage.filters import (
+    gaussian,
+    threshold_otsu,
+    median,
+    threshold_multiotsu,
+)
+
+from brainio import brainio
+from neuro.generic_neuro_tools import (
+    save_brain,
+    transform_all_channels_to_standard_space,
+)
+
+
+def get_fiber_track(
+    registered_atlas_path,
+    source_image_path,
+    seed_point,
+    output_path=None,
+    normalising_factor=4,
+    erosion_selem=(5, 5, 5),
+):
+
+    """
+    gets segmentation image of optical fibers using the background channel and a seed-point,
+    tested on 200um and 400um fibers
+
+    :param erosion_selem:
+    :param normalising_factor:
+    :param output_path:
+    :param registered_atlas_path:
+    :param source_image_path:
+    :param seed_point:
+    :return:
+    """
+    brain_mask = brainio.load_any(str(registered_atlas_path)) != 0
+    brain = brainio.load_any(str(source_image_path))
+
+    brain_median_filtered = median(brain)
+    otsu_threshold = threshold_multiotsu(brain_median_filtered)[0]
+    brain_segmentation = (
+        brain_median_filtered < otsu_threshold / normalising_factor
+    ) * brain_mask
+
+    segmentation_eroded = binary_erosion(
+        brain_segmentation, selem=erosion_selem
+    )
+    segmentation_eroded_dilated = binary_dilation(segmentation_eroded)
+    fiber_track_image = flood(
+        segmentation_eroded_dilated.astype(np.int16), seed_point
+    )
+
+    if output_path is not None:
+        save_brain(fiber_track_image, source_image_path, output_path)
+    else:
+        return fiber_track_image
+
+
+def get_lesion(
+    reg_dir,
+    allen_structure_id=672,
+    erosion_selem=np.ones((5, 5, 5)),
+    lesion_threshold=0.5,
+    sigma=5,
+    minimum_object_size=8000,
+    otsu_scale_factor=2,
+    transformed_brain_file_fname="registered_downsampled_channel_0.nii",
+    atlas_fname="annotations.nii",
+    output_fname="lesion_mask.nii",
+):
+    """
+
+    :param atlas_fname:
+    :param otsu_scale_factor:
+    :param transformed_brain_file_fname:
+    :param minimum_object_size: any blobs below this size will be ignored
+    :param sigma: sigma for gaussian filtering
+    :param lesion_threshold:
+    :param reg_dir: directory containing cellfinder/amap registration output files
+    :param allen_structure_id: the id of the structure that contains the lesion
+    :param erosion_selem:
+    :return:
+    """
+    reg_dir = pathlib.Path(reg_dir)
+    annotations_path = reg_dir / atlas_fname
+    source_image_path = reg_dir / transformed_brain_file_fname
+
+    if not os.path.isfile(str(source_image_path)):
+        print("channel not in std space... gotta transform")
+        transform_all_channels_to_standard_space(reg_dir)
+
+    brain_mask = brainio.load_any(str(annotations_path)) == allen_structure_id
+    brain = brainio.load_any(str(source_image_path))
+
+    brain_otsu_thresh = threshold_otsu(brain)
+    brain = brain < (brain_otsu_thresh / otsu_scale_factor)
+    brain = brain * brain_mask
+    brain = gaussian(brain, sigma) > lesion_threshold
+    brain = binary_erosion(brain, selem=erosion_selem)
+    brain = skimage.morphology.label(brain)
+    biggest = skimage.morphology.remove_small_objects(
+        brain, min_size=minimum_object_size
+    )
+    save_brain(biggest, annotations_path, reg_dir / output_fname)
+
+
+def get_lesion_sizes(reg_dir):
+    reg_dir = pathlib.Path(reg_dir)
+    labels = brainio.load_any(str(reg_dir / "lesion_mask.nii"))
+    all_structures_ids = np.unique(labels)
+    return [np.count_nonzero(labels == idx) for idx in all_structures_ids]
diff --git a/neuro/points/points_to_brainrender.py b/neuro/points/points_to_brainrender.py
index b2d0794..e8a6b55 100644
--- a/neuro/points/points_to_brainrender.py
+++ b/neuro/points/points_to_brainrender.py
@@ -8,6 +8,8 @@
 import argparse
 import imlib.IO.cells as cells_io
 from imlib.general.numerical import check_positive_float, check_positive_int
+from imlib.general.system import ensure_directory_exists
+from pathlib import Path
 
 
 def run(
@@ -34,6 +36,9 @@ def run(
 
     cells["x"] = max_z - cells["x"]
 
+    print("Ensuring output directory exists")
+    ensure_directory_exists(output_filename.parent)
+
     print(f"Saving to: {output_filename}")
     cells.to_hdf(output_filename, key=key, mode="w")
 
@@ -51,7 +56,7 @@ def get_parser():
     )
     parser.add_argument(
         dest="output_filename",
-        type=str,
+        type=Path,
         help="Output filename. Should end with '.h5'",
     )
 
diff --git a/neuro/visualise/__init__.py b/neuro/visualise/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/neuro/visualise/amap_vis.py b/neuro/visualise/amap_vis.py
new file mode 100644
index 0000000..565d915
--- /dev/null
+++ b/neuro/visualise/amap_vis.py
@@ -0,0 +1,126 @@
+import napari
+from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
+from pathlib import Path
+from imlib.IO.structures import load_structures_as_df
+from imlib.anatomy.structures.structures_tree import (
+    atlas_value_to_name,
+    UnknownAtlasValue,
+)
+from neuro.atlas_tools.paths import Paths
+from imlib.source.source_files import get_structures_path
+from neuro.visualise.vis_tools import (
+    display_raw,
+    display_downsampled,
+    display_registration,
+)
+
+
+def parser():
+    parser = ArgumentParser(formatter_class=ArgumentDefaultsHelpFormatter)
+    parser = cli_parse(parser)
+    return parser
+
+
+def cli_parse(parser):
+    cli_parser = parser.add_argument_group("Visualisation options")
+
+    cli_parser.add_argument(
+        dest="amap_directory",
+        type=str,
+        help="Path to the amap output directory..",
+    )
+
+    cli_parser.add_argument(
+        "-r" "--raw",
+        dest="raw",
+        action="store_true",
+        help="Display the raw image (as a virtual stack) "
+        "rather than the downsampled",
+    )
+    cli_parser.add_argument(
+        "-c",
+        "--raw-channels",
+        dest="raw_channels",
+        type=str,
+        nargs="+",
+        help="Paths to N additional raw channels to view. Will only work if "
+        "using the raw image viewer.",
+    )
+    cli_parser.add_argument(
+        "-m" "--memory",
+        dest="memory",
+        action="store_true",
+        help="Load data into RAM. ",
+    )
+
+    return parser
+
+
+def main():
+    print("Starting amap viewer")
+    args = parser().parse_args()
+
+    structures_path = get_structures_path()
+    structures_df = load_structures_as_df(structures_path)
+
+    if not args.memory:
+        print(
+            "By default amap_vis does not load data into memory. "
+            "To speed up visualisation, use the '-m' flag. Be aware "
+            "this will make the viewer slower to open initially."
+        )
+
+    paths = Paths(args.amap_directory)
+    with napari.gui_qt():
+        v = napari.Viewer(title="amap viewer")
+        if (
+            Path(paths.registered_atlas_path).exists()
+            and Path(paths.boundaries_file_path).exists()
+        ):
+
+            if args.raw:
+                image_scales = display_raw(v, args)
+
+            else:
+                if Path(paths.downsampled_brain_path).exists():
+                    image_scales = display_downsampled(v, args, paths)
+
+                else:
+                    raise FileNotFoundError(
+                        f"The downsampled image: "
+                        f"{paths.downsampled_brain_path} could not be found. "
+                        f"Please ensure this is the correct "
+                        f"directory and that amap has completed. "
+                    )
+
+            labels = display_registration(
+                v,
+                paths.registered_atlas_path,
+                paths.boundaries_file_path,
+                image_scales,
+                memory=args.memory,
+            )
+
+            @labels.mouse_move_callbacks.append
+            def get_connected_component_shape(layer, event):
+                val = layer.get_value()
+                if val != 0 and val is not None:
+                    try:
+                        region = atlas_value_to_name(val, structures_df)
+                        msg = f"{region}"
+                    except UnknownAtlasValue:
+                        msg = "Unknown region"
+                else:
+                    msg = "No label here!"
+                layer.help = msg
+
+        else:
+            raise FileNotFoundError(
+                f"The directory: '{args.amap_directory}' does not "
+                f"appear to be complete. Please ensure this is the correct "
+                f"directory and that amap has completed."
+            )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/neuro/visualise/vis_tools.py b/neuro/visualise/vis_tools.py
new file mode 100644
index 0000000..035e736
--- /dev/null
+++ b/neuro/visualise/vis_tools.py
@@ -0,0 +1,209 @@
+from pathlib import Path
+
+import numpy as np
+from brainio import brainio
+from natsort import natsorted
+
+from napari.utils.io import magic_imread
+
+from imlib.general.config import get_config_obj
+from imlib.general.system import get_sorted_file_paths, get_text_lines
+from vispy.color import Colormap
+
+label_red = Colormap([[0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0]])
+
+
+def prepare_load_nii(nii_path, memory=False):
+    """
+    Transforms a nii file into the same coordinate space as the raw data
+    :param nii_path: Path to the nii file
+    :param memory: Load data into memory
+    :return: Numpy array in the correct coordinate space
+    """
+    nii_path = str(nii_path)
+    image = brainio.load_any(nii_path, as_numpy=memory)
+    image = np.swapaxes(image, 2, 0)
+    return image
+
+
+def load_additional_downsampled_images(
+    viewer,
+    amap_directory,
+    paths,
+    search_string="downsampled_",
+    extension=".nii",
+    memory=False,
+):
+    """
+    Loads additional downsampled (i.e. from nii) images into napari viewer
+    :param viewer: Napari viewer object
+    :param amap_directory: Directory containing images
+    :param paths: amap paths object
+    :param search_string: String that defines the images.
+    Default: "downampled_"
+    :param extension: File extension of the downsampled images. Default: ".nii"
+    :param memory: Load data into memory
+    """
+
+    amap_directory = Path(amap_directory)
+
+    for file in amap_directory.iterdir():
+        if (
+            (file.suffix == ".nii")
+            and file.name.startswith(search_string)
+            and file != Path(paths.downsampled_brain_path)
+            and file != Path(paths.tmp__downsampled_filtered)
+        ):
+            print(
+                f"Found additional downsampled image: {file.name}, "
+                f"adding to viewer"
+            )
+            name = file.name.strip(search_string).strip(extension)
+            viewer.add_image(
+                prepare_load_nii(file, memory=memory), name=name,
+            )
+
+
+def add_raw_image(viewer, image_path, name):
+    """
+    Add a raw image (as a virtual stack) to the napari viewer
+    :param viewer: Napari viewer object
+    :param image_path: Path to the raw data
+    :param str name: Name to give the data
+    """
+    paths = get_sorted_file_paths(image_path, file_extension=".tif")
+    images = magic_imread(paths, use_dask=True, stack=True)
+    viewer.add_image(images, name=name, opacity=0.6, blending="additive")
+
+
+def get_image_scales(log_entries, config_file):
+    """
+    Returns the scaling from downsampled data to raw data
+    :param log_entries: Entries parsed from the log file
+    :param config_file: Path to the amap config file
+    :return: Tuple of scaling factors
+    """
+    config_obj = get_config_obj(config_file)
+    atlas_conf = config_obj["atlas"]
+    pixel_sizes = atlas_conf["pixel_size"]
+    x_scale = float(pixel_sizes["x"]) / float(log_entries["x_pixel_um"])
+    y_scale = float(pixel_sizes["y"]) / float(log_entries["y_pixel_um"])
+    z_scale = float(pixel_sizes["z"]) / float(log_entries["z_pixel_um"])
+    return z_scale, y_scale, x_scale
+
+
+def display_raw(viewer, args):
+    """
+    Display raw data
+    :param viewer:
+    :param args:
+    :return:
+    """
+    print(
+        "Starting raw image viewer. Streaming full-resolution data."
+        " This may be slow."
+    )
+    log_entries = read_log_file(get_most_recent_log(args.amap_directory))
+    config_file = Path(args.amap_directory, "config.conf")
+    image_scales = get_image_scales(log_entries, config_file)
+    add_raw_image(viewer, log_entries["image_paths"], name="Raw data")
+
+    if args.raw_channels:
+        for raw_image in args.raw_channels:
+            name = Path(raw_image).name
+            print(f"Found additional raw image to add to viewer: " f"{name}")
+            add_raw_image(viewer, raw_image, name=name)
+
+    return image_scales
+
+
+def display_downsampled(viewer, args, paths):
+    """
+    Display downsampled data
+    :param viewer:
+    :param args:
+    :param paths:
+    :return:
+    """
+    image_scales = (1, 1, 1)
+    load_additional_downsampled_images(
+        viewer, args.amap_directory, paths, memory=args.memory
+    )
+
+    viewer.add_image(
+        prepare_load_nii(paths.downsampled_brain_path, memory=args.memory),
+        name="Downsampled raw data",
+    )
+
+    return image_scales
+
+
+def display_registration(
+    viewer, atlas, boundaries, image_scales, memory=False
+):
+    """
+    Display results of the registration
+    :param viewer: napari viewer object
+    :param atlas: Annotations in sample space
+    :param boundaries: Annotation boundaries in sample space
+    :param tuple image_scales: Scaling of images from annotations -> data
+    :param memory: Load data into memory
+    """
+    viewer.add_image(
+        prepare_load_nii(boundaries, memory=memory),
+        name="Outlines",
+        contrast_limits=[0, 1],
+        colormap=("label_red", label_red),
+        scale=image_scales,
+    )
+
+    # labels added last so on top
+    labels = viewer.add_labels(
+        prepare_load_nii(atlas, memory=memory),
+        name="Annotations",
+        opacity=0.2,
+        scale=image_scales,
+    )
+    return labels
+
+
+def get_most_recent_log(directory, log_pattern="amap*.log"):
+    """
+    Returns the most recent amap log file (for parsing of arguments)
+    :param directory:
+    :param log_pattern: String pattern that defines the log
+    :return: Path to the most recent log file
+    """
+    directory = Path(directory)
+    return natsorted(directory.glob(log_pattern))[-1]
+
+
+def read_log_file(
+    log_file,
+    log_entries_to_get=[
+        "x_pixel_um",
+        "y_pixel_um",
+        "z_pixel_um",
+        "image_paths",
+        "registration_config",
+    ],
+    separator=": ",
+):
+    """
+    Reads an amap log file, and returns a dict of entries corresponding to
+    "log_entries_to_get"
+    :param log_file: Path to the log file
+    :param log_entries_to_get: List of strings corresponding to entries
+    in the log file
+    :param separator: Separator between the log item label and the entry.
+    Default: ": "
+    :return: A dict of the entries and labels
+    """
+    lines = get_text_lines(log_file)
+    entries = {}
+    for line in lines:
+        for entry in log_entries_to_get:
+            if line.startswith(entry):
+                entries[entry] = line.strip(entry + separator)
+
+    return entries
diff --git a/setup.py b/setup.py
index feee2ff..f290e0c 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 requirements = [
     "numpy",
     "scikit-image",
-    "pandas",
+    "pandas<=0.25.3,>=0.25.1",
     "napari",
     "brainrender",
     "imlib",
@@ -49,6 +49,7 @@
             "points_to_brainrender = "
             "neuro.points.points_to_brainrender:main",
             "heatmap = neuro.heatmap.heatmap:cli",
+            "amap_vis = neuro.visualise.amap_vis:main",
         ]
     },
     zip_safe=False,
diff --git a/tests/data/heatmap/heatmap.nii b/tests/data/heatmap/heatmap.nii
index 72ef829..2694417 100644
Binary files a/tests/data/heatmap/heatmap.nii and b/tests/data/heatmap/heatmap.nii differ