diff --git a/.dockerignore b/.dockerignore
index b5bfff59..a51f661d 100644
--- a/.dockerignore
+++ b/.dockerignore
@@ -8,8 +8,8 @@ build/**/*
 build
 dist/**/*
 dist
-fmriprep.egg-info/**/*
-fmriprep.egg-info
+dmriprep.egg-info/**/*
+dmriprep.egg-info
 .eggs/**/*
 .eggs
 
@@ -39,4 +39,4 @@ out/
 .zenodo.json
 .travis.yml
 .readthedocs.yml
-CONTRIBUTING.md
\ No newline at end of file
+CONTRIBUTING.rst
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 00000000..cd5a1171
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1 @@
+dmriprep/_version.py export-subst
diff --git a/.github/config.yml b/.github/config.yml
index c65ed7d8..70c14c71 100644
--- a/.github/config.yml
+++ b/.github/config.yml
@@ -5,8 +5,8 @@ newPRWelcomeComment: >
   Thanks for opening this pull request!
   We have detected this is the first time for you to contribute
   to *dMRIPrep*.
-  Please check out our [contributing guidelines](https://github.com/nipreps/dmriprep/blob/master/CONTRIBUTING.md).
-  
+  Please check out our [contributing guidelines](https://github.com/nipreps/dmriprep/blob/master/CONTRIBUTING.rst).
+
   We invite you to list yourself as a *dMRIPrep* contributor, so if your name
   is not already mentioned, please modify the
   [``.zenodo.json``](https://github.com/nipreps/dmriprep/blob/master/.zenodo.json)
@@ -25,4 +25,4 @@ newPRWelcomeComment: >
   Of course, if you want to opt-out this time there is no
   problem at all with adding your name later.
   You will be always welcome to add it in the future whenever
-  you feel it should be listed.
\ No newline at end of file
+  you feel it should be listed.
diff --git a/.github/release-drafter.yml b/.github/release-drafter.yml
new file mode 100644
index 00000000..075f95b6
--- /dev/null
+++ b/.github/release-drafter.yml
@@ -0,0 +1,5 @@
+template: |
+  ## Release Notes
+
+  ## CHANGES
+  $CHANGES
diff --git a/.github/stale.yml b/.github/stale.yml
new file mode 100644
index 00000000..28e0d5c1
--- /dev/null
+++ b/.github/stale.yml
@@ -0,0 +1,20 @@
+# Number of days of inactivity before an issue becomes stale
+daysUntilStale: 900
+# Number of days of inactivity before a stale issue is closed
+daysUntilClose: 200
+# Issues with these labels will never be considered stale
+exemptLabels:
+  - pinned
+  - security
+  - feature
+  - help wanted
+  - low priority
+# Label to use when marking an issue as stale
+staleLabel: stale
+# Comment to post when marking an issue as stale. Set to `false` to disable
+markComment: >
+  This issue has been automatically marked as stale because it has not had
+  recent activity. It will be closed if no further activity occurs. Thank you
+  for your contributions.
+# Comment to post when closing a stale issue. Set to `false` to disable
+closeComment: false
diff --git a/.gitignore b/.gitignore
index a8237ba3..3cae9a46 100644
--- a/.gitignore
+++ b/.gitignore
@@ -103,6 +103,3 @@ ENV/
 
 # Mac OS nonsense:
 .DS_Store
-
-#kubernetes stuff
-kubernetes/jobs/
diff --git a/.readthedocs.yml b/.readthedocs.yml
new file mode 100644
index 00000000..14d721c1
--- /dev/null
+++ b/.readthedocs.yml
@@ -0,0 +1,20 @@
+# .readthedocs.yml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+  configuration: docs/conf.py
+
+# Optionally build your docs in additional formats such as PDF and ePub
+formats: all
+
+python:
+  install:
+    - method: pip
+      path: .
+      extra_requirements:
+        - doc
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 00000000..293b1af1
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,35 @@
+# Config file for automatic testing at travis-ci.org
+
+language: python
+matrix:
+  include:
+    - python: 3.5
+      dist: trusty
+      sudo: false
+    - python: 3.6
+      dist: trusty
+      sudo: false
+    - python: 3.7
+      dist: xenial
+      sudo: true
+
+# Command to install dependencies, e.g. pip install -r requirements.txt --use-mirrors
+install: pip install -U tox-travis
+
+# Command to run tests, e.g. python setup.py test
+script: tox
+
+# Assuming you have installed the travis-ci CLI tool, after you
+# create the Github repo and add it to Travis, run the
+# following command to finish PyPI deployment setup:
+# $ travis encrypt --add deploy.password
+deploy:
+  provider: pypi
+  distributions: sdist bdist_wheel
+  user: nipy
+  password:
+    secure: PLEASE_REPLACE_ME
+  on:
+    tags: true
+    repo: tigrlab/dmriprep
+    python: 3.7
diff --git a/.zenodo.json b/.zenodo.json
index f823aa44..615c6db0 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -39,6 +39,17 @@
     {
       "name": "nrajamani3",
       "type": "Researcher"
+    },
+    {
+        "name": "Mansour, Salim",
+        "affiliation": "The Centre for Addiction and Mental Health",
+        "type": "Researcher"
+    },
+    {
+        "name": "Forde, Natalie",
+        "affiliation": "The Centre for Addiction and Mental Health",
+        "orcid": "0000-0002-2005-8774",
+        "type": "Researcher"
     }
   ],
   "keywords": [
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
new file mode 100644
index 00000000..956ccfc1
--- /dev/null
+++ b/CONTRIBUTING.rst
@@ -0,0 +1,181 @@
+.. highlight:: shell
+
+============
+Contributing
+============
+
+We love contributions! dmriprep is open source, built on open source,
+and we'd love to have you hang out in our community.
+
+**Imposter syndrome disclaimer**: We want your help. No, really.
+
+There may be a little voice inside your head that is telling you that
+you're not ready to be an open source contributor; that your skills
+aren't nearly good enough to contribute. What could you possibly offer a
+project like this one?
+
+We assure you - the little voice in your head is wrong. If you can
+write code at all, you can contribute code to open source. Contributing
+to open source projects is a fantastic way to advance one's coding
+skills. Writing perfect code isn't the measure of a good developer (that
+would disqualify all of us!); it's trying to create something, making
+mistakes, and learning from those mistakes. That's how we all improve,
+and we are happy to help others learn.
+
+Being an open source contributor doesn't just mean writing code, either.
+You can help out by writing documentation, tests, or even giving
+feedback about the project (and yes - that includes giving feedback
+about the contribution process). Some of these contributions may be the
+most valuable to the project as a whole, because you're coming to the
+project with fresh eyes, so you can see the errors and assumptions that
+seasoned contributors have glossed over.
+
+Installing a development version of dmriprep
+--------------------------------------------
+
+First, you can install a development version of dmriprep by cloning this repository
+and then typing::
+
+    $ pip install -e .[dev]
+
+Activate the pre-commit formatting hook by typing::
+
+    $ pre-commit install
+
+Before committing your work, you can check for formatting issues or error by typing::
+
+    $ make lint
+    $ make test
+
+Types of Contributions
+----------------------
+
+You can contribute in many ways:
+
+Report Bugs
+~~~~~~~~~~~
+
+Report bugs at https://github.com/nipreps/dmriprep/issues.
+
+If you are reporting a bug, please include:
+
+* Your operating system name and version.
+* Any details about your local setup that might be helpful in troubleshooting.
+* Detailed steps to reproduce the bug.
+
+Fix Bugs
+~~~~~~~~
+
+Look through the GitHub issues for bugs. Anything tagged with "bug" and "help
+wanted" is open to whoever wants to implement it.
+
+Implement Features
+~~~~~~~~~~~~~~~~~~
+
+Look through the GitHub issues for features. Anything tagged with "enhancement"
+and "help wanted" is open to whoever wants to implement it.
+
+Write Documentation
+~~~~~~~~~~~~~~~~~~~
+
+dmriprep could always use more documentation, whether as part of the
+official dmriprep docs, in docstrings, or even on the web in blog posts,
+articles, and such.
+
+Submit Feedback
+~~~~~~~~~~~~~~~
+
+The best way to send feedback is to file an issue at https://github.com/nipreps/dmriprep/issues.
+
+If you are proposing a feature:
+
+* Explain in detail how it would work.
+* Keep the scope as narrow as possible, to make it easier to implement.
+* Remember that this is a volunteer-driven project, and that contributions
+  are welcome :)
+
+Get Started!
+------------
+
+Ready to contribute? Here's how to set up `dmriprep` for local development.
+
+1. Fork the `dmriprep` repo on GitHub.
+2. Clone your fork locally::
+
+    $ git clone git@github.com:your_name_here/dmriprep.git
+
+3. Install your local copy into a virtualenv. Assuming you have virtualenvwrapper installed, this is how you set up your fork for local development::
+
+    $ mkvirtualenv dmriprep
+    $ cd dmriprep/
+    $ python setup.py develop
+
+4. Create a branch for local development::
+
+    $ git checkout -b name-of-your-bugfix-or-feature
+
+   Now you can make your changes locally.
+
+5. When you're done making changes, check that your changes pass flake8 and the
+   tests, including testing other Python versions with tox::
+
+    $ flake8 dmriprep tests
+    $ python setup.py test or py.test
+    $ tox
+
+   To get flake8 and tox, just pip install them into your virtualenv.
+
+6. Commit your changes and push your branch to GitHub::
+
+    $ git add .
+    $ git commit -m "Your detailed description of your changes."
+    $ git push origin name-of-your-bugfix-or-feature
+
+7. Submit a pull request through the GitHub website.
+
+Pull Request Guidelines
+-----------------------
+
+Before you submit a pull request, check that it meets these guidelines:
+
+1. The pull request should include tests.
+2. If the pull request adds functionality, the docs should be updated. Put
+   your new functionality into a function with a docstring, and add the
+   feature to the list in README.rst.
+3. The pull request should work for Python 3.5, 3.6 and 3.7, and for PyPy. Check
+   https://travis-ci.org/tigrlab/dmriprep/pull_requests
+   and make sure that the tests pass for all supported Python versions.
+
+When opening a pull request, please use one of the following prefixes:
+
+* **[ENH]** for enhancements
+* **[FIX]** for bug fixes
+* **[TST]** for new or updated tests
+* **[DOC]** for new or updated documentation
+* **[STY]** for stylistic changes
+* **[REF]** for refactoring existing code
+
+Tips
+----
+
+To run a subset of tests::
+
+$ py.test tests.test_dmriprep
+
+
+Deploying
+---------
+
+A reminder for the maintainers on how to deploy.
+Make sure all your changes are committed (including an entry in HISTORY.rst).
+Then run::
+
+$ bumpversion patch # possible: major / minor / patch
+$ git push
+$ git push --tags
+
+Travis will then deploy to PyPI if tests pass.
+
+The imposter syndrome disclaimer was originally written by
+`Adrienne Lowe <https://github.com/adriennefriend>`_ for a `PyCon talk <https://www.youtube.com/watch?v=6Uj746j9Heo>`_, and was
+adapted based on its use in the README file for the `MetPy project <https://github.com/Unidata/MetPy>`_.
diff --git a/Dockerfile b/Dockerfile
index 9c707115..fe16af09 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -170,7 +170,7 @@ RUN pip install --no-cache-dir "templateflow>=0.4.0,<0.5.0a0" && \
     find $HOME/.cache/templateflow -type d -exec chmod go=u {} + && \
     find $HOME/.cache/templateflow -type f -exec chmod go=u {} +
 
-# Installing FMRIPREP
+# Installing DMRIPREP
 COPY . /src/dmriprep
 ARG VERSION
 # Force static versioning within container
diff --git a/Makefile b/Makefile
new file mode 100644
index 00000000..03b6b514
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,94 @@
+.PHONY: clean clean-test clean-pyc clean-build docs help
+.DEFAULT_GOAL := help
+
+define BROWSER_PYSCRIPT
+import os, webbrowser, sys
+
+try:
+	from urllib import pathname2url
+except:
+	from urllib.request import pathname2url
+
+webbrowser.open("file://" + pathname2url(os.path.abspath(sys.argv[1])))
+endef
+export BROWSER_PYSCRIPT
+
+define PRINT_HELP_PYSCRIPT
+import re, sys
+
+for line in sys.stdin:
+	match = re.match(r'^([a-zA-Z_-]+):.*?## (.*)$$', line)
+	if match:
+		target, help = match.groups()
+		print("%-20s %s" % (target, help))
+endef
+export PRINT_HELP_PYSCRIPT
+
+BROWSER := python -c "$$BROWSER_PYSCRIPT"
+
+help:
+	@python -c "$$PRINT_HELP_PYSCRIPT" < $(MAKEFILE_LIST)
+
+clean: clean-build clean-pyc clean-test ## remove all build, test, coverage and Python artifacts
+
+clean-build: ## remove build artifacts
+	rm -fr build/
+	rm -fr dist/
+	rm -fr .eggs/
+	find . -name '*.egg-info' -exec rm -fr {} +
+	find . -name '*.egg' -exec rm -f {} +
+
+clean-pyc: ## remove Python file artifacts
+	find . -name '*.pyc' -exec rm -f {} +
+	find . -name '*.pyo' -exec rm -f {} +
+	find . -name '*~' -exec rm -f {} +
+	find . -name '__pycache__' -exec rm -fr {} +
+
+clean-test: ## remove test and coverage artifacts
+	rm -fr .tox/
+	rm -f .coverage
+	rm -fr htmlcov/
+	rm -fr .pytest_cache
+
+lint: ## check style with flake8
+	flake8 dmriprep tests
+
+test: ## run tests quickly with the default Python
+	py.test
+
+test-all: ## run tests on every Python version with tox
+	tox
+
+coverage: ## check code coverage quickly with the default Python
+	coverage run --source dmriprep -m pytest
+	coverage report -m
+	coverage html
+	$(BROWSER) htmlcov/index.html
+
+docs: ## generate Sphinx HTML documentation, including API docs
+	rm -f docs/dmriprep*.rst
+	rm -f docs/modules.rst
+	sphinx-apidoc -o docs/ dmriprep
+	$(MAKE) -C docs clean
+	$(MAKE) -C docs html
+	$(BROWSER) docs/_build/html/index.html
+
+servedocs: docs ## compile the docs watching for changes
+	watchmedo shell-command -p '*.rst' -c '$(MAKE) -C docs html' -R -D .
+
+release: dist ## package and upload a release
+	twine upload dist/*
+
+dist: clean ## builds source and wheel package
+	python setup.py sdist
+	python setup.py bdist_wheel
+	ls -l dist
+
+install: clean ## install the package to the active Python's site-packages
+	python setup.py install
+
+docker: docker
+	docker build --rm -t nipreps/dmriprep:latest \
+		--build-arg BUILD_DATE=`date -u +"%Y-%m-%dT%H:%M:%SZ"` \
+		--build-arg VCS_REF=`git rev-parse --short HEAD` \
+		--build-arg VERSION=$( python get_version.py ) .
diff --git a/README.rst b/README.rst
index f4add394..3070c2fd 100644
--- a/README.rst
+++ b/README.rst
@@ -4,13 +4,10 @@ dmriprep
 
 .. image:: https://badgen.net/badge/chat/on%20mattermost/blue
    :target: https://mattermost.brainhack.org/brainhack/channels/dmriprep
-   
+
 .. image:: https://img.shields.io/pypi/v/dmriprep.svg
         :target: https://pypi.python.org/pypi/dmriprep
 
-.. image:: https://circleci.com/gh/nipreps/dmriprep.svg?style=svg
-    :target: https://circleci.com/gh/nipreps/dmriprep
-
 .. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3392201.svg
    :target: https://doi.org/10.5281/zenodo.3392201
 
@@ -49,4 +46,4 @@ segmentation, skullstripping etc.) providing outputs that can be
 easily submitted to a variety of tractography algorithms.
 
 [Documentation `dmriprep.org <https://dmriprep.readthedocs.io>`_]
-[Support `neurostars.org <https://neurostars.org/tags/fmriprep>`_]
+[Support `neurostars.org <https://neurostars.org/tags/dmriprep>`_]
diff --git a/dmriprep/__init__.py b/dmriprep/__init__.py
index 72e61e5a..35bac792 100644
--- a/dmriprep/__init__.py
+++ b/dmriprep/__init__.py
@@ -15,7 +15,7 @@
     '__packagename__',
 ]
 
-# cmp is not used by fmriprep, so ignore nipype-generated warnings
+# cmp is not used by dmriprep, so ignore nipype-generated warnings
 _warnings.filterwarnings('ignore', r'cmp not installed')
 _warnings.filterwarnings('ignore', r'This has not been fully tested. Please report any failures.')
 _warnings.filterwarnings('ignore', r"can't resolve package from __spec__ or __package__")
diff --git a/dmriprep/cli/run.py b/dmriprep/cli/run.py
index 05f74b0c..30cb3246 100755
--- a/dmriprep/cli/run.py
+++ b/dmriprep/cli/run.py
@@ -75,6 +75,8 @@ def get_parser():
     # Re-enable when option is actually implemented
     # g_bids.add_argument('-r', '--run-id', action='store', default='single_run',
     #                     help='select a specific run to be processed')
+    g_bids.add_argument('-a', '--acq-id', action='store',
+                        help='select a specific acquisition to be processed')
 
     g_perfm = parser.add_argument_group('Options to handle performance')
     g_perfm.add_argument('--nprocs', '--n_cpus', '-n-cpus', action='store', type=int,
@@ -97,7 +99,8 @@ def get_parser():
 
     g_conf = parser.add_argument_group('Workflow configuration')
     g_conf.add_argument(
-        '--ignore', required=False, action='store', nargs="+", default=[], choices=['sdc'],
+        '--ignore', required=False, action='store', nargs="+", default=[],
+        choices=['denoising', 'unringing', 'fieldmaps'],
         help='ignore selected aspects of the input dataset to disable corresponding '
              'parts of the workflow (a space delimited list)')
     g_conf.add_argument(
@@ -129,6 +132,17 @@ def get_parser():
                         help='do not use a random seed for skull-stripping - will ensure '
                              'run-to-run replicability when used with --omp-nthreads 1')
 
+    # Dwi options
+    g_dwi = parser.add_argument_group('Specific options for handling dwi scans')
+    g_dwi.add_argument('--concat-dwis', required=False, action='store',
+                       nargs="+", default=[],
+                       help='space-delimited list of acq-<label>')
+    g_dwi.add_argument('--b0-thresh', required=False, action='store', default=5,
+                       help='Threshold for b0 value')
+    g_dwi.add_argument('--output-resolution', required=False, action='store',
+                       help='isotropic voxel size in mm the data will be resampled '
+                            'to before running eddy')
+
     # Fieldmap options
     g_fmap = parser.add_argument_group('Specific options for handling fieldmaps')
     g_fmap.add_argument('--fmap-bspline', action='store_true', default=False,
@@ -151,6 +165,15 @@ def get_parser():
         help='Path to FreeSurfer license key file. Get it (for free) by registering'
              ' at https://surfer.nmr.mgh.harvard.edu/registration.html')
 
+    # Eddy options
+    g_eddy = parser.add_argument_group('Specific options for FSL eddy preprocessing')
+    g_eddy.add_argument('--eddy-config', action='store', type=Path,
+                        help='Path to json file with input arguments for running '
+                             'eddy')
+    g_eddy.add_argument('--acqp-file', action='store', type=Path,
+                        help='Path to acquisition parameters file for topup/eddy '
+                             'instead of generating it from the dwi json')
+
     # Surface generation xor
     g_surfs = parser.add_argument_group('Surface preprocessing options')
     g_surfs_xor = g_surfs.add_mutually_exclusive_group()
@@ -546,8 +569,14 @@ def build_workflow(opts, retval):
     )
 
     retval['workflow'] = init_dmriprep_wf(
+        acq_id=opts.acq_id,
+        acqp_file=opts.acqp_file,
         anat_only=opts.anat_only,
+        b0_thresh=opts.b0_thresh,
+        concat_dwis=opts.concat_dwis,
         debug=opts.debug,
+        eddy_config=opts.eddy_config,
+        fmap_bspline=opts.fmap_bspline,
         force_syn=opts.force_syn,
         freesurfer=opts.run_reconall,
         hires=opts.hires,
@@ -557,6 +586,7 @@ def build_workflow(opts, retval):
         low_mem=opts.low_mem,
         omp_nthreads=omp_nthreads,
         output_dir=str(output_dir),
+        output_resolution=opts.output_resolution,
         output_spaces=output_spaces,
         run_uuid=run_uuid,
         skull_strip_fixed_seed=opts.skull_strip_fixed_seed,
diff --git a/dmriprep/data/eddy_parameters.json b/dmriprep/data/eddy_parameters.json
new file mode 100644
index 00000000..0c7b08f9
--- /dev/null
+++ b/dmriprep/data/eddy_parameters.json
@@ -0,0 +1,27 @@
+{
+    "flm": "quadratic",
+    "slm": "none",
+    "fwhm": 0,
+    "niter": 5,
+    "fep": false,
+    "interp": "spline",
+    "method": "jac",
+    "nvoxhp": 1000,
+    "fudge_factor": 10,
+    "dont_sep_offs_move": false,
+    "dont_peas": false,
+    "repol": true,
+    "outlier_nstd": 4,
+    "outlier_nvox": 250,
+    "outlier_type": "sw",
+    "outlier_pos": false,
+    "outlier_sqr": false,
+    "estimate_mov_by_susceptibility": false,
+    "is_shelled": false,
+    "num_threads": 1,
+    "use_cuda": false,
+    "cnr_maps": true,
+    "residuals": true,
+    "output_type": "NIFTI_GZ",
+    "args": ""
+}
diff --git a/dmriprep/due.py b/dmriprep/due.py
new file mode 100644
index 00000000..9a1c4dd0
--- /dev/null
+++ b/dmriprep/due.py
@@ -0,0 +1,74 @@
+# emacs: at the end of the file
+# ex: set sts=4 ts=4 sw=4 et:
+# ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #
+"""
+
+Stub file for a guaranteed safe import of duecredit constructs:  if duecredit
+is not available.
+
+To use it, place it into your project codebase to be imported, e.g. copy as
+
+    cp stub.py /path/tomodule/module/due.py
+
+Note that it might be better to avoid naming it duecredit.py to avoid shadowing
+installed duecredit.
+
+Then use in your code as
+
+    from .due import due, Doi, BibTeX, Text
+
+See  https://github.com/duecredit/duecredit/blob/master/README.md for examples.
+
+Origin:     Originally a part of the duecredit
+Copyright:  2015-2019  DueCredit developers
+License:    BSD-2
+"""
+
+__version__ = '0.0.8'
+
+
+class InactiveDueCreditCollector(object):
+    """Just a stub at the Collector which would not do anything"""
+    def _donothing(self, *args, **kwargs):
+        """Perform no good and no bad"""
+        pass
+
+    def dcite(self, *args, **kwargs):
+        """If I could cite I would"""
+        def nondecorating_decorator(func):
+            return func
+        return nondecorating_decorator
+
+    active = False
+    activate = add = cite = dump = load = _donothing
+
+    def __repr__(self):
+        return self.__class__.__name__ + '()'
+
+
+def _donothing_func(*args, **kwargs):
+    """Perform no good and no bad"""
+    pass
+
+
+try:
+    from duecredit import due, BibTeX, Doi, Url, Text
+    if 'due' in locals() and not hasattr(due, 'cite'):
+        raise RuntimeError(
+            "Imported due lacks .cite. DueCredit is now disabled")
+except Exception as e:
+    if not isinstance(e, ImportError):
+        import logging
+        logging.getLogger("duecredit").error(
+            "Failed to import duecredit due to %s" % str(e))
+    # Initiate due stub
+    due = InactiveDueCreditCollector()
+    BibTeX = Doi = Url = Text = _donothing_func
+
+# Emacs mode definitions
+# Local Variables:
+# mode: python
+# py-indent-offset: 4
+# tab-width: 4
+# indent-tabs-mode: nil
+# End:
diff --git a/dmriprep/interfaces/fmap.py b/dmriprep/interfaces/fmap.py
new file mode 100644
index 00000000..85c0996b
--- /dev/null
+++ b/dmriprep/interfaces/fmap.py
@@ -0,0 +1,612 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+Interfaces to deal with the various types of fieldmap sources
+"""
+
+import numpy as np
+import nibabel as nb
+from nipype import logging
+from nipype.utils.filemanip import fname_presuffix
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
+    SimpleInterface, InputMultiObject)
+
+LOGGER = logging.getLogger('nipype.interface')
+
+
+class FieldEnhanceInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input fieldmap')
+    in_mask = File(exists=True, desc='brain mask')
+    in_magnitude = File(exists=True, desc='input magnitude')
+    unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap')
+    despike = traits.Bool(True, usedefault=True, desc='run despike filter')
+    bspline_smooth = traits.Bool(True, usedefault=True, desc='run 3D bspline smoother')
+    mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations')
+    despike_threshold = traits.Float(0.2, usedefault=True, desc='mask erosion iterations')
+    num_threads = traits.Int(1, usedefault=True, nohash=True, desc='number of jobs')
+
+
+class FieldEnhanceOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    out_unwrapped = File(desc='unwrapped fieldmap')
+
+
+class FieldEnhance(SimpleInterface):
+    """
+    The FieldEnhance interface wraps a workflow to massage the input fieldmap
+    and return it masked, despiked, etc.
+    """
+    input_spec = FieldEnhanceInputSpec
+    output_spec = FieldEnhanceOutputSpec
+
+    def _run_interface(self, runtime):
+        from scipy import ndimage as sim
+
+        fmap_nii = nb.load(self.inputs.in_file)
+        data = np.squeeze(fmap_nii.get_data().astype(np.float32))
+
+        # Despike / denoise (no-mask)
+        if self.inputs.despike:
+            data = _despike2d(data, self.inputs.despike_threshold)
+
+        mask = None
+        if isdefined(self.inputs.in_mask):
+            masknii = nb.load(self.inputs.in_mask)
+            mask = masknii.get_data().astype(np.uint8)
+
+            # Dilate mask
+            if self.inputs.mask_erode > 0:
+                struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1)
+                mask = sim.binary_erosion(
+                    mask, struc,
+                    iterations=self.inputs.mask_erode
+                ).astype(np.uint8)  # pylint: disable=no-member
+
+        self._results['out_file'] = fname_presuffix(
+            self.inputs.in_file, suffix='_enh', newpath=runtime.cwd)
+        datanii = nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header)
+
+        if self.inputs.unwrap:
+            data = _unwrap(data, self.inputs.in_magnitude, mask)
+            self._results['out_unwrapped'] = fname_presuffix(
+                self.inputs.in_file, suffix='_unwrap', newpath=runtime.cwd)
+            nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header).to_filename(
+                self._results['out_unwrapped'])
+
+        if not self.inputs.bspline_smooth:
+            datanii.to_filename(self._results['out_file'])
+            return runtime
+        else:
+            from ..utils import bspline as fbsp
+            from statsmodels.robust.scale import mad
+
+            # Fit BSplines (coarse)
+            bspobj = fbsp.BSplineFieldmap(datanii, weights=mask,
+                                          njobs=self.inputs.num_threads)
+            bspobj.fit()
+            smoothed1 = bspobj.get_smoothed()
+
+            # Manipulate the difference map
+            diffmap = data - smoothed1.get_data()
+            sderror = mad(diffmap[mask > 0])
+            LOGGER.info('SD of error after B-Spline fitting is %f', sderror)
+            errormask = np.zeros_like(diffmap)
+            errormask[np.abs(diffmap) > (10 * sderror)] = 1
+            errormask *= mask
+
+            nslices = 0
+            try:
+                errorslice = np.squeeze(np.argwhere(errormask.sum(0).sum(0) > 0))
+                nslices = errorslice[-1] - errorslice[0]
+            except IndexError:  # mask is empty, do not refine
+                pass
+
+            if nslices > 1:
+                diffmapmsk = mask[..., errorslice[0]:errorslice[-1]]
+                diffmapnii = nb.Nifti1Image(
+                    diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk,
+                    datanii.affine, datanii.header)
+
+                bspobj2 = fbsp.BSplineFieldmap(diffmapnii, knots_zooms=[24., 24., 4.],
+                                               njobs=self.inputs.num_threads)
+                bspobj2.fit()
+                smoothed2 = bspobj2.get_smoothed().get_data()
+
+                final = smoothed1.get_data().copy()
+                final[..., errorslice[0]:errorslice[-1]] += smoothed2
+            else:
+                final = smoothed1.get_data()
+
+            nb.Nifti1Image(final, datanii.affine, datanii.header).to_filename(
+                self._results['out_file'])
+
+        return runtime
+
+
+class FieldToRadSInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input fieldmap')
+    fmap_range = traits.Float(desc='range of input field map')
+
+
+class FieldToRadSOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    fmap_range = traits.Float(desc='range of input field map')
+
+
+class FieldToRadS(SimpleInterface):
+    """
+    The FieldToRadS converts from arbitrary units to rad/s
+    """
+    input_spec = FieldToRadSInputSpec
+    output_spec = FieldToRadSOutputSpec
+
+    def _run_interface(self, runtime):
+        fmap_range = None
+        if isdefined(self.inputs.fmap_range):
+            fmap_range = self.inputs.fmap_range
+        self._results['out_file'], self._results['fmap_range'] = _torads(
+            self.inputs.in_file, fmap_range, newpath=runtime.cwd)
+        return runtime
+
+
+class FieldToHzInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input fieldmap')
+    range_hz = traits.Float(mandatory=True, desc='range of input field map')
+
+
+class FieldToHzOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+
+
+class FieldToHz(SimpleInterface):
+    """
+    The FieldToHz converts from arbitrary units to Hz
+    """
+    input_spec = FieldToHzInputSpec
+    output_spec = FieldToHzOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = _tohz(
+            self.inputs.in_file, self.inputs.range_hz, newpath=runtime.cwd)
+        return runtime
+
+
+class Phasediff2FieldmapInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input fieldmap')
+    metadata = traits.Dict(mandatory=True, desc='BIDS metadata dictionary')
+
+
+class Phasediff2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+
+
+class Phasediff2Fieldmap(SimpleInterface):
+    """
+    Convert a phase difference map into a fieldmap in Hz
+    """
+    input_spec = Phasediff2FieldmapInputSpec
+    output_spec = Phasediff2FieldmapOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = phdiff2fmap(
+            self.inputs.in_file,
+            _delta_te(self.inputs.metadata),
+            newpath=runtime.cwd)
+        return runtime
+
+
+class Phases2FieldmapInputSpec(BaseInterfaceInputSpec):
+    phase_files = InputMultiObject(
+        File(exists=True), mandatory=True, desc='list of phase1, phase2 files')
+    metadatas = traits.List(
+        traits.Dict, mandatory=True, desc='list of phase1, phase2 metadata dicts')
+
+
+class Phases2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    phasediff_metadata = traits.Dict(desc='the phasediff metadata')
+
+
+class Phases2Fieldmap(SimpleInterface):
+    """
+    Convert a phase1, phase2 into a difference map
+    """
+    input_spec = Phases2FieldmapInputSpec
+    output_spec = Phases2FieldmapOutputSpec
+
+    def _run_interface(self, runtime):
+        # Get the echo times
+        fmap_file, merged_metadata = phases2fmap(self.inputs.phase_files, self.inputs.metadatas,
+                                                 newpath=runtime.cwd)
+        self._results['phasediff_metadata'] = merged_metadata
+        self._results['out_file'] = fmap_file
+        return runtime
+
+
+def _despike2d(data, thres, neigh=None):
+    """
+    despiking as done in FSL fugue
+    """
+
+    if neigh is None:
+        neigh = [-1, 0, 1]
+    nslices = data.shape[-1]
+
+    for k in range(nslices):
+        data2d = data[..., k]
+
+        for i in range(data2d.shape[0]):
+            for j in range(data2d.shape[1]):
+                vals = []
+                thisval = data2d[i, j]
+                for ii in neigh:
+                    for jj in neigh:
+                        try:
+                            vals.append(data2d[i + ii, j + jj])
+                        except IndexError:
+                            pass
+                vals = np.array(vals)
+                patch_range = vals.max() - vals.min()
+                patch_med = np.median(vals)
+
+                if (patch_range > 1e-6 and
+                        (abs(thisval - patch_med) / patch_range) > thres):
+                    data[i, j, k] = patch_med
+    return data
+
+
+def _unwrap(fmap_data, mag_file, mask=None):
+    from math import pi
+    from nipype.interfaces.fsl import PRELUDE
+    magnii = nb.load(mag_file)
+
+    if mask is None:
+        mask = np.ones_like(fmap_data, dtype=np.uint8)
+
+    fmapmax = max(abs(fmap_data[mask > 0].min()), fmap_data[mask > 0].max())
+    fmap_data *= pi / fmapmax
+
+    nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz')
+    nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz')
+    nb.Nifti1Image(magnii.get_data(), magnii.affine).to_filename('fmap_mag.nii.gz')
+
+    # Run prelude
+    res = PRELUDE(phase_file='fmap_rad.nii.gz',
+                  magnitude_file='fmap_mag.nii.gz',
+                  mask_file='fmap_mask.nii.gz').run()
+
+    unwrapped = nb.load(res.outputs.unwrapped_phase_file).get_data() * (fmapmax / pi)
+    return unwrapped
+
+
+def get_ees(in_meta, in_file=None):
+    """
+    Calculate the *effective echo spacing* :math:`t_\\text{ees}`
+    for an input :abbr:`EPI (echo-planar imaging)` scan.
+
+
+    There are several procedures to calculate the effective
+    echo spacing. The basic one is that an ``EffectiveEchoSpacing``
+    field is set in the JSON sidecar. The following examples
+    use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the
+    j-axis encoding direction.
+
+    >>> meta = {'EffectiveEchoSpacing': 0.00059,
+    ...         'PhaseEncodingDirection': 'j-'}
+    >>> get_ees(meta)
+    0.00059
+
+    If the *total readout time* :math:`T_\\text{ro}` (``TotalReadoutTime``
+    BIDS field) is provided, then the effective echo spacing can be
+    calculated reading the number of voxels :math:`N_\\text{PE}` along the
+    readout direction and the parallel acceleration
+    factor of the EPI
+
+      .. math ::
+
+           =  T_\\text{ro} \\,  (N_\\text{PE} / f_\\text{acc} - 1)^{-1}
+
+    where :math:`N_y` is the number of pixels along the phase-encoding direction
+    :math:`y`, and :math:`f_\\text{acc}` is the parallel imaging acceleration factor
+    (:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`,
+    :abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.).
+
+    >>> meta = {'TotalReadoutTime': 0.02596,
+    ...         'PhaseEncodingDirection': 'j-',
+    ...         'ParallelReductionFactorInPlane': 2}
+    >>> get_ees(meta, in_file='epi.nii.gz')
+    0.00059
+
+    Some vendors, like Philips, store different parameter names
+    (see http://dbic.dartmouth.edu/pipermail/mrusers/attachments/\
+20141112/eb1d20e6/attachment.pdf):
+
+    >>> meta = {'WaterFatShift': 8.129,
+    ...         'MagneticFieldStrength': 3,
+    ...         'PhaseEncodingDirection': 'j-',
+    ...         'ParallelReductionFactorInPlane': 2}
+    >>> get_ees(meta, in_file='epi.nii.gz')
+    0.00041602630141921826
+
+    """
+
+    import nibabel as nb
+    from fmriprep.interfaces.fmap import _get_pe_index
+
+    # Use case 1: EES is defined
+    ees = in_meta.get('EffectiveEchoSpacing', None)
+    if ees is not None:
+        return ees
+
+    # All other cases require the parallel acc and npe (N vox in PE dir)
+    acc = float(in_meta.get('ParallelReductionFactorInPlane', 1.0))
+    npe = nb.load(in_file).shape[_get_pe_index(in_meta)]
+    etl = npe // acc
+
+    # Use case 2: TRT is defined
+    trt = in_meta.get('TotalReadoutTime', None)
+    if trt is not None:
+        return trt / (etl - 1)
+
+    # Use case 3 (philips scans)
+    wfs = in_meta.get('WaterFatShift', None)
+    if wfs is not None:
+        fstrength = in_meta['MagneticFieldStrength']
+        wfd_ppm = 3.4  # water-fat diff in ppm
+        g_ratio_mhz_t = 42.57  # gyromagnetic ratio for proton (1H) in MHz/T
+        wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
+        return wfs / (wfs_hz * etl)
+
+    raise ValueError('Unknown effective echo-spacing specification')
+
+
+def get_trt(in_meta, in_file=None):
+    """
+    Calculate the *total readout time* for an input
+    :abbr:`EPI (echo-planar imaging)` scan.
+
+
+    There are several procedures to calculate the total
+    readout time. The basic one is that a ``TotalReadoutTime``
+    field is set in the JSON sidecar. The following examples
+    use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the
+    j-axis encoding direction.
+
+    >>> meta = {'TotalReadoutTime': 0.02596}
+    >>> get_trt(meta)
+    0.02596
+
+    If the *effective echo spacing* :math:`t_\\text{ees}`
+    (``EffectiveEchoSpacing`` BIDS field) is provided, then the
+    total readout time can be calculated reading the number
+    of voxels along the readout direction :math:`T_\\text{ro}`
+    and the parallel acceleration factor of the EPI :math:`f_\\text{acc}`.
+
+      .. math ::
+
+          T_\\text{ro} = t_\\text{ees} \\, (N_\\text{PE} / f_\\text{acc} - 1)
+
+    >>> meta = {'EffectiveEchoSpacing': 0.00059,
+    ...         'PhaseEncodingDirection': 'j-',
+    ...         'ParallelReductionFactorInPlane': 2}
+    >>> get_trt(meta, in_file='epi.nii.gz')
+    0.02596
+
+    Some vendors, like Philips, store different parameter names:
+
+    >>> meta = {'WaterFatShift': 8.129,
+    ...         'MagneticFieldStrength': 3,
+    ...         'PhaseEncodingDirection': 'j-',
+    ...         'ParallelReductionFactorInPlane': 2}
+    >>> get_trt(meta, in_file='epi.nii.gz')
+    0.018721183563864822
+
+    """
+
+    # Use case 1: TRT is defined
+    trt = in_meta.get('TotalReadoutTime', None)
+    if trt is not None:
+        return trt
+
+    # All other cases require the parallel acc and npe (N vox in PE dir)
+    acc = float(in_meta.get('ParallelReductionFactorInPlane', 1.0))
+    npe = nb.load(in_file).shape[_get_pe_index(in_meta)]
+    etl = npe // acc
+
+    # Use case 2: TRT is defined
+    ees = in_meta.get('EffectiveEchoSpacing', None)
+    if ees is not None:
+        return ees * (etl - 1)
+
+    # Use case 3 (philips scans)
+    wfs = in_meta.get('WaterFatShift', None)
+    if wfs is not None:
+        fstrength = in_meta['MagneticFieldStrength']
+        wfd_ppm = 3.4  # water-fat diff in ppm
+        g_ratio_mhz_t = 42.57  # gyromagnetic ratio for proton (1H) in MHz/T
+        wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
+        return wfs / wfs_hz
+
+    raise ValueError('Unknown total-readout time specification')
+
+
+def _get_pe_index(meta):
+    pe = meta['PhaseEncodingDirection']
+    try:
+        return {'i': 0, 'j': 1, 'k': 2}[pe[0]]
+    except KeyError:
+        raise RuntimeError('"%s" is an invalid PE string' % pe)
+
+
+def _torads(in_file, fmap_range=None, newpath=None):
+    """
+    Convert a field map to rad/s units
+
+    If fmap_range is None, the range of the fieldmap
+    will be automatically calculated.
+
+    Use fmap_range=0.5 to convert from Hz to rad/s
+    """
+    from math import pi
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    out_file = fname_presuffix(in_file, suffix='_rad', newpath=newpath)
+    fmapnii = nb.load(in_file)
+    fmapdata = fmapnii.get_data()
+
+    if fmap_range is None:
+        fmap_range = max(abs(fmapdata.min()), fmapdata.max())
+    fmapdata = fmapdata * (pi / fmap_range)
+    out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header)
+    out_img.set_data_dtype('float32')
+    out_img.to_filename(out_file)
+    return out_file, fmap_range
+
+
+def _tohz(in_file, range_hz, newpath=None):
+    """Convert a field map to Hz units"""
+    from math import pi
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    out_file = fname_presuffix(in_file, suffix='_hz', newpath=newpath)
+    fmapnii = nb.load(in_file)
+    fmapdata = fmapnii.get_data()
+    fmapdata = fmapdata * (range_hz / pi)
+    out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header)
+    out_img.set_data_dtype('float32')
+    out_img.to_filename(out_file)
+    return out_file
+
+
+def phdiff2fmap(in_file, delta_te, newpath=None):
+    r"""
+    Converts the input phase-difference map into a fieldmap in Hz,
+    using the eq. (1) of [Hutton2002]_:
+
+    .. math::
+
+        \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}}
+
+
+    In this case, we do not take into account the gyromagnetic ratio of the
+    proton (:math:`\gamma`), since it will be applied inside TOPUP:
+
+    .. math::
+
+        \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}}
+
+    """
+    import math
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+    #  GYROMAG_RATIO_H_PROTON_MHZ = 42.576
+
+    out_file = fname_presuffix(in_file, suffix='_fmap', newpath=newpath)
+    image = nb.load(in_file)
+    data = (image.get_data().astype(np.float32) / (2. * math.pi * delta_te))
+    nii = nb.Nifti1Image(data, image.affine, image.header)
+    nii.set_data_dtype(np.float32)
+    nii.to_filename(out_file)
+    return out_file
+
+
+def phases2fmap(phase_files, metadatas, newpath=None):
+    """Calculates a phasediff from two phase images. Assumes monopolar
+    readout. """
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+    from copy import deepcopy
+
+    phasediff_file = fname_presuffix(phase_files[0], suffix='_phasediff', newpath=newpath)
+    echo_times = [meta.get("EchoTime") for meta in metadatas]
+    if None in echo_times or echo_times[0] == echo_times[1]:
+        raise RuntimeError()
+    # Determine the order of subtraction
+    short_echo_index = echo_times.index(min(echo_times))
+    long_echo_index = echo_times.index(max(echo_times))
+
+    short_phase_image = phase_files[short_echo_index]
+    long_phase_image = phase_files[long_echo_index]
+
+    image0 = nb.load(short_phase_image)
+    phase0 = image0.get_fdata()
+    image1 = nb.load(long_phase_image)
+    phase1 = image1.get_fdata()
+
+    def rescale_image(img):
+        if np.any(img < -128):
+            # This happens sometimes on 7T fieldmaps
+            LOGGER.info("Found negative values in phase image: rescaling")
+            imax = img.max()
+            imin = img.min()
+            scaled = 2 * ((img - imin) / (imax - imin) - 0.5)
+            return np.pi * scaled
+        mask = img > 0
+        imax = img.max()
+        imin = img.min()
+        max_check = imax - 4096
+        if np.abs(max_check) > 10 or np.abs(imin) > 10:
+            LOGGER.warning("Phase image may be scaled incorrectly: check results")
+        return mask * (img / 2048 * np.pi - np.pi)
+
+    # Calculate fieldmaps
+    rad0 = rescale_image(phase0)
+    rad1 = rescale_image(phase1)
+    a = np.cos(rad0)
+    b = np.sin(rad0)
+    c = np.cos(rad1)
+    d = np.sin(rad1)
+    fmap = -np.arctan2(b * c - a * d, a * c + b * d)
+
+    phasediff_nii = nb.Nifti1Image(fmap, image0.affine)
+    phasediff_nii.set_data_dtype(np.float32)
+    phasediff_nii.to_filename(phasediff_file)
+
+    merged_metadata = deepcopy(metadatas[0])
+    del merged_metadata['EchoTime']
+    merged_metadata['EchoTime1'] = float(echo_times[short_echo_index])
+    merged_metadata['EchoTime2'] = float(echo_times[long_echo_index])
+
+    return phasediff_file, merged_metadata
+
+
+def _delta_te(in_values, te1=None, te2=None):
+    r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict"""
+    if isinstance(in_values, float):
+        te2 = in_values
+        te1 = 0.
+
+    if isinstance(in_values, dict):
+        te1 = in_values.get('EchoTime1')
+        te2 = in_values.get('EchoTime2')
+
+        if not all((te1, te2)):
+            te2 = in_values.get('EchoTimeDifference')
+            te1 = 0
+
+    if isinstance(in_values, list):
+        te2, te1 = in_values
+        if isinstance(te1, list):
+            te1 = te1[1]
+        if isinstance(te2, list):
+            te2 = te2[1]
+
+    # For convienience if both are missing we should give one error about them
+    if te1 is None and te2 is None:
+        raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found. '
+                           'Please consult the BIDS specification.')
+    if te1 is None:
+        raise RuntimeError(
+            'EchoTime1 metadata field not found. Please consult the BIDS specification.')
+    if te2 is None:
+        raise RuntimeError(
+            'EchoTime2 metadata field not found. Please consult the BIDS specification.')
+
+    return abs(float(te2) - float(te1))
diff --git a/dmriprep/interfaces/fsl.py b/dmriprep/interfaces/fsl.py
new file mode 100644
index 00000000..5b78c6bb
--- /dev/null
+++ b/dmriprep/interfaces/fsl.py
@@ -0,0 +1,467 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+Interface for FSL Eddy (temporary until updated Eddy spec is merged)
+"""
+from __future__ import print_function, division, unicode_literals, \
+    absolute_import
+from builtins import str
+
+import os
+
+from nipype.utils.filemanip import fname_presuffix
+
+from nipype.interfaces.base import (traits, TraitedSpec, File, isdefined)
+from nipype.interfaces.fsl.base import FSLCommand, FSLCommandInputSpec
+
+
+class EddyInputSpec(FSLCommandInputSpec):
+    in_file = File(
+        exists=True,
+        mandatory=True,
+        argstr='--imain=%s',
+        desc='File containing all the images to estimate distortions for')
+    in_mask = File(
+        exists=True,
+        mandatory=True,
+        argstr='--mask=%s',
+        desc='Mask to indicate brain')
+    in_index = File(
+        exists=True,
+        mandatory=True,
+        argstr='--index=%s',
+        desc=('File containing indices for all volumes in --imain '
+              'into --acqp and --topup'))
+    in_acqp = File(
+        exists=True,
+        mandatory=True,
+        argstr='--acqp=%s',
+        desc='File containing acquisition parameters')
+    in_bvec = File(
+        exists=True,
+        mandatory=True,
+        argstr='--bvecs=%s',
+        desc='File containing the b-vectors for all volumes in --imain')
+    in_bval = File(
+        exists=True,
+        mandatory=True,
+        argstr='--bvals=%s',
+        desc='File containing the b-values for all volumes in --imain')
+    out_base = traits.Str(
+        default_value='eddy_corrected',
+        usedefault=True,
+        argstr='--out=%s',
+        desc='basename for output (warped) image')
+    session = File(
+        exists=True,
+        argstr='--session=%s',
+        desc='File containing session indices for all volumes in --imain')
+    in_topup_fieldcoef = File(
+        exists=True,
+        argstr='--topup=%s',
+        requires=['in_topup_movpar'],
+        desc='topup file containing the field coefficients')
+    in_topup_movpar = File(
+        exists=True,
+        requires=['in_topup_fieldcoef'],
+        desc='topup movpar.txt file')
+    field = File(
+        exists=True,
+        argstr='--field=%s',
+        desc=('NonTOPUP fieldmap scaled in Hz - filename has to be provided '
+              'without an extension. TOPUP is strongly recommended'))
+    field_mat = File(
+        exists=True,
+        argstr='--field_mat=%s',
+        requires=['field'],
+        desc=('Matrix that specifies the relative locations of the field '
+              'specified by --field and first volume in file --imain'))
+    flm = traits.Enum(
+        'linear',
+        'quadratic',
+        'cubic',
+        argstr='--flm=%s',
+        desc='First level EC model')
+    slm = traits.Enum(
+        'none',
+        'linear',
+        'quadratic',
+        argstr='--slm=%s',
+        desc='Second level EC model')
+    fwhm = traits.Float(
+        default_value=0,
+        usedefault=True,
+        argstr='--fwhm=%s',
+        desc='FWHM for conditioning filter when estimating the parameters')
+    niter = traits.Int(
+        default_value=5,
+        usedefault=True,
+        argstr='--niter=%s',
+        desc='Number of iterations')
+    fep = traits.Bool(
+        False, argstr='--fep', desc='Fill empty planes in x- or y-directions')
+    initrand = traits.Bool(
+        False,
+        argstr='--initrand',
+        desc='Resets rand for when selecting voxels',
+        min_ver='5.0.10')
+    interp = traits.Enum(
+        'spline',
+        'trilinear',
+        argstr='--interp=%s',
+        desc='Interpolation model for estimation step')
+    method = traits.Enum(
+        'jac',
+        'lsr',
+        argstr='--resamp=%s',
+        desc='Final resampling method (jacobian/least squares)')
+    nvoxhp = traits.Int(
+        default_value=1000,
+        usedefault=True,
+        argstr='--nvoxhp=%s',
+        desc='# of voxels used to estimate the hyperparameters')
+    fudge_factor = traits.Float(
+        default_value=10.0,
+        usedefault=True,
+        argstr='--ff=%s',
+        desc='Fudge factor for hyperparameter error variance')
+    dont_sep_offs_move = traits.Bool(
+        False,
+        argstr='--dont_sep_offs_move',
+        desc='Do NOT attempt to separate field offset from subject movement')
+    dont_peas = traits.Bool(
+        False,
+        argstr='--dont_peas',
+        desc='Do NOT perform a post-eddy alignment of shells')
+    repol = traits.Bool(
+        default_value=False,
+        argstr='--repol',
+        desc='Detect and replace outlier slices')
+    outlier_nstd = traits.Int(
+        argstr='--ol_nstd',
+        desc='Number of std off to qualify as outlier',
+        requires=['repol'],
+        min_ver='5.0.10')
+    outlier_nvox = traits.Int(
+        argstr='--ol_nvox',
+        desc='Min # of voxels in a slice for inclusion in outlier detection',
+        requires=['repol'],
+        min_ver='5.0.10')
+    outlier_type = traits.Enum(
+        'sw',
+        'gw',
+        'both',
+        argstr='--ol_type',
+        desc='Type of outliers, slicewise (sw), groupwise (gw) or both (both)',
+        requires=['repol'],
+        min_ver='5.0.10')
+    outlier_pos = traits.Bool(
+        False,
+        argstr='--ol_pos',
+        desc='Consider both positive and negative outliers if set',
+        requires=['repol'],
+        min_ver='5.0.10')
+    outlier_sqr = traits.Bool(
+        False,
+        argstr='--ol_sqr',
+        desc='Consider outliers among sums-of-squared differences if set',
+        requires=['repol'],
+        min_ver='5.0.10')
+    mb = traits.Int(
+        argstr='--mb=%s',
+        desc='Multi-band factor',
+        min_ver='5.0.10')
+    mb_offs = traits.Enum(
+        0,
+        1,
+        -1,
+        argstr='--mb_offs=%s',
+        desc=('Multi-band offset (-1 if bottom slice removed, 1 if '
+              'top slice removed'),
+        requires=['mb'],
+        min_ver='5.0.10')
+    mporder = traits.Int(
+        argstr='--mporder=%s',
+        desc='Order of slice-to-vol movement model',
+        requires=['slspec'],
+        min_ver='5.0.11')
+    s2v_niter = traits.Int(
+        argstr='--s2v_niter=%s',
+        desc='Number of iterations for slice-to-vol',
+        requires=['slspec'],
+        min_ver='5.0.11')
+    s2v_lambda = traits.Int(
+        agstr='--s2v_lambda',
+        desc='Regularisation weight for slice-to-vol movement (reasonable range 1-10)',
+        requires=['slspec'],
+        min_ver='5.0.11')
+    s2v_interp = traits.Enum(
+        'trilinear',
+        'spline',
+        argstr='--s2v_interp=%s',
+        desc='Slice-to-vol interpolation model for estimation step',
+        requires=['slspec'],
+        min_ver='5.0.11')
+    slspec = traits.File(
+        exists=True,
+        argstr='--slspec=%s',
+        desc='Name of text file completely specifying slice/group acquisition',
+        xor=['json'],
+        min_ver='5.0.11')
+    json = traits.File(
+        exists=True,
+        argstr='--json=%s',
+        desc='Name of .json text file with information about slice timing',
+        xor=['slspec'],
+        min_ver='6.0.1')
+    estimate_move_by_susceptibility = traits.Bool(
+        False,
+        argstr='--estimate_move_by_susceptibility',
+        desc='Estimate how susceptibility field changes with subject movement',
+        min_ver='6.0.1')
+    mbs_niter = traits.Int(
+        argstr='--mbs_niter=%s',
+        desc='Number of iterations for MBS estimation',
+        requires=['estimate_move_by_susceptibility'],
+        min_ver='6.0.1')
+    mbs_lambda = traits.Int(
+        argstr='--mbs_lambda=%s',
+        desc='Weighting of regularisation for MBS estimation',
+        requires=['estimate_move_by_susceptibility'],
+        min_ver='6.0.1')
+    mbs_ksp = traits.Int(
+        argstr='--mbs_ksp=%smm',
+        desc='Knot-spacing for MBS field estimation',
+        requires=['estimate_move_by_susceptibility'],
+        min_ver='6.0.1')
+    is_shelled = traits.Bool(
+        False,
+        argstr='--data_is_shelled',
+        desc=('Override internal check to ensure that data are acquired '
+              'on a set of b-value shells'))
+    num_threads = traits.Int(
+        default_value=1,
+        usedefault=True,
+        nohash=True,
+        desc="Number of openmp threads to use")
+    use_cuda = traits.Bool(False, desc="Run eddy using cuda gpu")
+    cnr_maps = traits.Bool(
+        False, desc='Output CNR-Maps', argstr='--cnr_maps', min_ver='5.0.10')
+    residuals = traits.Bool(
+        False, desc='Output Residuals', argstr='--residuals', min_ver='5.0.10')
+
+
+class EddyOutputSpec(TraitedSpec):
+    out_corrected = File(
+        exists=True, desc='4D image file containing all the corrected volumes')
+    out_parameter = File(
+        exists=True,
+        desc=('Text file with parameters defining the field and'
+              'movement for each scan'))
+    out_rotated_bvecs = File(
+        exists=True, desc='File containing rotated b-values for all volumes')
+    out_movement_rms = File(
+        exists=True, desc='Summary of the "total movement" in each volume')
+    out_restricted_movement_rms = File(
+        exists=True,
+        desc=('Summary of the "total movement" in each volume '
+              'disregarding translation in the PE direction'))
+    out_shell_alignment_parameters = File(
+        exists=True,
+        desc=('Text file containing rigid body movement parameters '
+              'between the different shells as estimated by a '
+              'post-hoc mutual information based registration'))
+    out_shell_pe_translation_parameters = File(
+        exists=True,
+        desc=('Text file containing translation along the PE-direction '
+              'between the different shells as estimated by a '
+              'post-hoc mutual information based registration'))
+    out_outlier_map = File(
+        exists=True,
+        desc=('Matrix where rows represent volumes and columns represent '
+              'slices. "0" indicates that scan-slice is not an outlier '
+              'and "1" indicates that it is'))
+    out_outlier_n_stdev_map = File(
+        exists=True,
+        desc=('Matrix where rows represent volumes and columns represent '
+              'slices. Values indicate number of standard deviations off the '
+              'mean difference between observation and prediction is'))
+    out_outlier_n_sqr_stdev_map = File(
+        exists=True,
+        desc=('Matrix where rows represent volumes and columns represent '
+              'slices. Values indicate number of standard deivations off the '
+              'square root of the mean squared difference between observation '
+              'and prediction is'))
+    out_outlier_report = File(
+        exists=True,
+        desc=('Text file with a plain language report on what '
+              'outlier slices eddy has found'))
+    out_outlier_free = File(
+        exists=True,
+        desc=('4D image file not corrected for susceptibility or eddy-'
+              'current distortions or subject movement but with outlier '
+              'slices replaced'))
+    out_movement_over_time = File(
+        exists=True,
+        desc=('Text file containing translations (mm) and rotations '
+              '(radians) for each excitation'))
+    out_cnr_maps = File(
+        exists=True, desc='path/name of file with the cnr_maps')
+    out_residuals = File(
+        exists=True, desc='path/name of file with the residuals')
+
+
+class Eddy(FSLCommand):
+    """
+    Interface for FSL eddy, a tool for estimating and correcting eddy
+    currents induced distortions. `User guide
+    <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide>`_ and
+    `more info regarding acqp file
+    <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/Faq#How_do_I_know_what_to_put_into_my_--acqp_file.3F>`_.
+    Examples
+    --------
+    >>> from nipype.interfaces.fsl import Eddy
+    >>> eddy = Eddy()
+    >>> eddy.inputs.in_file = 'epi.nii'
+    >>> eddy.inputs.in_mask  = 'epi_mask.nii'
+    >>> eddy.inputs.in_index = 'epi_index.txt'
+    >>> eddy.inputs.in_acqp  = 'epi_acqp.txt'
+    >>> eddy.inputs.in_bvec  = 'bvecs.scheme'
+    >>> eddy.inputs.in_bval  = 'bvals.scheme'
+    >>> eddy.inputs.use_cuda = True
+    >>> eddy.cmdline # doctest: +ELLIPSIS
+    'eddy_cuda --ff=10.0 --fwhm=0 --acqp=epi_acqp.txt --bvals=bvals.scheme \
+--bvecs=bvecs.scheme --imain=epi.nii --index=epi_index.txt \
+--mask=epi_mask.nii --niter=5 --nvoxhp=1000 --out=.../eddy_corrected'
+    >>> eddy.inputs.use_cuda = False
+    >>> eddy.cmdline # doctest: +ELLIPSIS
+    'eddy_openmp --ff=10.0 --fwhm=0 --acqp=epi_acqp.txt --bvals=bvals.scheme \
+--bvecs=bvecs.scheme --imain=epi.nii --index=epi_index.txt \
+--mask=epi_mask.nii --niter=5 --nvoxhp=1000 --out=.../eddy_corrected'
+    >>> res = eddy.run() # doctest: +SKIP
+    """
+    _cmd = 'eddy_openmp'
+    input_spec = EddyInputSpec
+    output_spec = EddyOutputSpec
+
+    _num_threads = 1
+
+    def __init__(self, **inputs):
+        super(Eddy, self).__init__(**inputs)
+        self.inputs.on_trait_change(self._num_threads_update, 'num_threads')
+        if not isdefined(self.inputs.num_threads):
+            self.inputs.num_threads = self._num_threads
+        else:
+            self._num_threads_update()
+        self.inputs.on_trait_change(self._use_cuda, 'use_cuda')
+        if isdefined(self.inputs.use_cuda):
+            self._use_cuda()
+
+    def _num_threads_update(self):
+        self._num_threads = self.inputs.num_threads
+        if not isdefined(self.inputs.num_threads):
+            if 'OMP_NUM_THREADS' in self.inputs.environ:
+                del self.inputs.environ['OMP_NUM_THREADS']
+        else:
+            self.inputs.environ['OMP_NUM_THREADS'] = str(
+                self.inputs.num_threads)
+
+    def _use_cuda(self):
+        self._cmd = 'eddy_cuda' if self.inputs.use_cuda else 'eddy_openmp'
+
+    def _run_interface(self, runtime):
+        # If 'eddy_openmp' is missing, use 'eddy'
+        FSLDIR = os.getenv('FSLDIR', '')
+        cmd = self._cmd
+        if all((FSLDIR != '', cmd == 'eddy_openmp',
+                not os.path.exists(os.path.join(FSLDIR, 'bin', cmd)))):
+            self._cmd = 'eddy'
+        runtime = super(Eddy, self)._run_interface(runtime)
+
+        # Restore command to avoid side-effects
+        self._cmd = cmd
+        return runtime
+
+    def _format_arg(self, name, spec, value):
+        if name == 'in_topup_fieldcoef':
+            return spec.argstr % value.split('_fieldcoef')[0]
+        if name == 'field':
+            return spec.argstr % fname_presuffix(value, use_ext=False)
+        if name == 'out_base':
+            return spec.argstr % os.path.abspath(value)
+        return super(Eddy, self)._format_arg(name, spec, value)
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['out_corrected'] = os.path.abspath(
+            '%s.nii.gz' % self.inputs.out_base)
+        outputs['out_parameter'] = os.path.abspath(
+            '%s.eddy_parameters' % self.inputs.out_base)
+
+        # File generation might depend on the version of EDDY
+        out_rotated_bvecs = os.path.abspath(
+            '%s.eddy_rotated_bvecs' % self.inputs.out_base)
+        out_movement_rms = os.path.abspath(
+            '%s.eddy_movement_rms' % self.inputs.out_base)
+        out_restricted_movement_rms = os.path.abspath(
+            '%s.eddy_restricted_movement_rms' % self.inputs.out_base)
+        out_shell_alignment_parameters = os.path.abspath(
+            '%s.eddy_post_eddy_shell_alignment_parameters' %
+            self.inputs.out_base)
+        out_shell_pe_translation_parameters = os.path.abspath(
+            '%s.eddy_post_eddy_shell_PE_translation_parameters' %
+            self.inputs.out_base)
+        out_outlier_map = os.path.abspath(
+            '%s.eddy_outlier_map' % self.inputs.out_base)
+        out_outlier_n_stdev_map = os.path.abspath(
+            '%s.eddy_outlier_n_stdev_map' % self.inputs.out_base)
+        out_outlier_n_sqr_stdev_map = os.path.abspath(
+            '%s.eddy_outlier_n_sqr_stdev_map' % self.inputs.out_base)
+        out_outlier_report = os.path.abspath(
+            '%s.eddy_outlier_report' % self.inputs.out_base)
+        if isdefined(self.inputs.repol) and self.inputs.repol:
+            out_outlier_free = os.path.abspath(
+                '%s.eddy_outlier_free_data' % self.inputs.out_base)
+            if os.path.exists(out_outlier_free):
+                outputs['out_outlier_free'] = out_outlier_free
+        if isdefined(self.inputs.mporder) and self.inputs.mporder > 0:
+            out_movement_over_time = os.path.abspath(
+                '%s.eddy_movement_over_time' % self.inputs.out_base)
+            if os.path.exists(out_movement_over_time):
+                outputs['out_movement_over_time'] = out_movement_over_time
+        if isdefined(self.inputs.cnr_maps) and self.inputs.cnr_maps:
+            out_cnr_maps = os.path.abspath(
+                '%s.eddy_cnr_maps.nii.gz' % self.inputs.out_base)
+            if os.path.exists(out_cnr_maps):
+                outputs['out_cnr_maps'] = out_cnr_maps
+        if isdefined(self.inputs.residuals) and self.inputs.residuals:
+            out_residuals = os.path.abspath(
+                '%s.eddy_residuals.nii.gz' % self.inputs.out_base)
+            if os.path.exists(out_residuals):
+                outputs['out_residuals'] = out_residuals
+
+        if os.path.exists(out_rotated_bvecs):
+            outputs['out_rotated_bvecs'] = out_rotated_bvecs
+        if os.path.exists(out_movement_rms):
+            outputs['out_movement_rms'] = out_movement_rms
+        if os.path.exists(out_restricted_movement_rms):
+            outputs['out_restricted_movement_rms'] = \
+                out_restricted_movement_rms
+        if os.path.exists(out_shell_alignment_parameters):
+            outputs['out_shell_alignment_parameters'] = \
+                out_shell_alignment_parameters
+        if os.path.exists(out_shell_pe_translation_parameters):
+            outputs['out_shell_pe_translation_parameters'] = \
+                out_shell_pe_translation_parameters
+        if os.path.exists(out_outlier_map):
+            outputs['out_outlier_map'] = out_outlier_map
+        if os.path.exists(out_outlier_n_stdev_map):
+            outputs['out_outlier_n_stdev_map'] = out_outlier_n_stdev_map
+        if os.path.exists(out_outlier_n_sqr_stdev_map):
+            outputs['out_outlier_n_sqr_stdev_map'] = out_outlier_n_sqr_stdev_map
+        if os.path.exists(out_outlier_report):
+            outputs['out_outlier_report'] = out_outlier_report
+
+        return outputs
diff --git a/dmriprep/interfaces/mrtrix3.py b/dmriprep/interfaces/mrtrix3.py
new file mode 100644
index 00000000..15c54628
--- /dev/null
+++ b/dmriprep/interfaces/mrtrix3.py
@@ -0,0 +1,254 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+Interfaces from Mrtrix3 (temporary until MRResize is merged)
+"""
+
+from nipype.interfaces.base import (traits, TraitedSpec, File)
+from nipype.interfaces.mrtrix3.base import MRTrix3BaseInputSpec, MRTrix3Base
+
+
+class DWIDenoiseInputSpec(MRTrix3BaseInputSpec):
+    in_file = File(
+        exists=True,
+        argstr='%s',
+        position=-2,
+        mandatory=True,
+        desc='input DWI image')
+    mask = File(
+        exists=True,
+        argstr='-mask %s',
+        position=1,
+        desc='mask image')
+    extent = traits.Tuple(
+        (traits.Int, traits.Int, traits.Int),
+        argstr='-extent %d,%d,%d',
+        desc='set the window size of the denoising filter. (default = 5,5,5)')
+    noise = File(
+        argstr='-noise %s',
+        name_template='%s_noise',
+        name_source='in_file',
+        keep_extension=True,
+        desc='the output noise map')
+    out_file = File(
+        argstr='%s',
+        position=-1,
+        name_template='%s_denoised',
+        name_source='in_file',
+        keep_extension=True,
+        desc='the output denoised DWI image')
+
+
+class DWIDenoiseOutputSpec(TraitedSpec):
+    noise = File(desc='the output noise map', exists=True)
+    out_file = File(desc='the output denoised DWI image', exists=True)
+
+
+class DWIDenoise(MRTrix3Base):
+    """
+    Denoise DWI data and estimate the noise level based on the optimal
+    threshold for PCA.
+
+    DWI data denoising and noise map estimation by exploiting data redundancy
+    in the PCA domain using the prior knowledge that the eigenspectrum of
+    random covariance matrices is described by the universal Marchenko Pastur
+    distribution.
+
+    Important note: image denoising must be performed as the first step of the
+    image processing pipeline. The routine will fail if interpolation or
+    smoothing has been applied to the data prior to denoising.
+
+    Note that this function does not correct for non-Gaussian noise biases.
+
+    For more information, see
+    <https://mrtrix.readthedocs.io/en/latest/reference/commands/dwidenoise.html>
+
+    Example
+    -------
+
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> denoise = mrt.DWIDenoise()
+    >>> denoise.inputs.in_file = 'dwi.mif'
+    >>> denoise.inputs.mask = 'mask.mif'
+    >>> denoise.inputs.noise = 'noise.mif'
+    >>> denoise.cmdline                               # doctest: +ELLIPSIS
+    'dwidenoise -mask mask.mif -noise noise.mif dwi.mif dwi_denoised.mif'
+    >>> denoise.run()                                 # doctest: +SKIP
+    """
+
+    _cmd = 'dwidenoise'
+    input_spec = DWIDenoiseInputSpec
+    output_spec = DWIDenoiseOutputSpec
+
+
+class MRDeGibbsInputSpec(MRTrix3BaseInputSpec):
+    in_file = File(
+        exists=True,
+        argstr='%s',
+        position=-2,
+        mandatory=True,
+        desc='input DWI image')
+    axes = traits.ListInt(
+        default_value=[0, 1],
+        usedefault=True,
+        sep=',',
+        minlen=2,
+        maxlen=2,
+        argstr='-axes %s',
+        desc='indicate the plane in which the data was acquired (axial = 0,1; '
+             'coronal = 0,2; sagittal = 1,2')
+    nshifts = traits.Int(
+        default_value=20,
+        usedefault=True,
+        argstr='-nshifts %d',
+        desc='discretization of subpixel spacing (default = 20)')
+    minW = traits.Int(
+        default_value=1,
+        usedefault=True,
+        argstr='-minW %d',
+        desc='left border of window used for total variation (TV) computation '
+             '(default = 1)')
+    maxW = traits.Int(
+        default_value=3,
+        usedefault=True,
+        argstr='-maxW %d',
+        desc='right border of window used for total variation (TV) computation '
+             '(default = 3)')
+    out_file = File(
+        name_template='%s_unr',
+        name_source='in_file',
+        keep_extension=True,
+        argstr='%s',
+        position=-1,
+        desc='the output unringed DWI image')
+
+
+class MRDeGibbsOutputSpec(TraitedSpec):
+    out_file = File(desc='the output unringed DWI image', exists=True)
+
+
+class MRDeGibbs(MRTrix3Base):
+    """
+    Remove Gibbs ringing artifacts.
+
+    This application attempts to remove Gibbs ringing artefacts from MRI images
+    using the method of local subvoxel-shifts proposed by Kellner et al.
+
+    This command is designed to run on data directly after it has been
+    reconstructed by the scanner, before any interpolation of any kind has
+    taken place. You should not run this command after any form of motion
+    correction (e.g. not after dwipreproc). Similarly, if you intend running
+    dwidenoise, you should run this command afterwards, since it has the
+    potential to alter the noise structure, which would impact on dwidenoise's
+    performance.
+
+    Note that this method is designed to work on images acquired with full
+    k-space coverage. Running this method on partial Fourier ('half-scan') data
+    may lead to suboptimal and/or biased results, as noted in the original
+    reference below. There is currently no means of dealing with this; users
+    should exercise caution when using this method on partial Fourier data, and
+    inspect its output for any obvious artefacts.
+
+    For more information, see
+    <https://mrtrix.readthedocs.io/en/latest/reference/commands/mrdegibbs.html>
+
+    Example
+    -------
+
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> unring = mrt.MRDeGibbs()
+    >>> unring.inputs.in_file = 'dwi.mif'
+    >>> unring.cmdline
+    'mrdegibbs -axes 0,1 -maxW 3 -minW 1 -nshifts 20 dwi.mif dwi_unr.mif'
+    >>> unring.run()                                 # doctest: +SKIP
+    """
+
+    _cmd = 'mrdegibbs'
+    input_spec = MRDeGibbsInputSpec
+    output_spec = MRDeGibbsOutputSpec
+
+
+class MRResizeInputSpec(MRTrix3BaseInputSpec):
+    in_file = File(
+        exists=True,
+        argstr='%s',
+        position=-2,
+        mandatory=True,
+        desc='input DWI image'
+    )
+    image_size = traits.Tuple(
+        (traits.Int, traits.Int, traits.Int),
+        argstr='-size %d,%d,%d',
+        mandatory=True,
+        desc='define the new image size for the output image. This should be '
+        'specified as a comma-separated list.',
+        xor=['voxel_size', 'scale_factor'],
+    )
+    voxel_size = traits.Tuple(
+        (traits.Float, traits.Float, traits.Float),
+        argstr='-voxel %d,%d,%d',
+        mandatory=True,
+        desc='define the new voxel size for the output image. This can be '
+        'specified either as a single value to be used for all '
+        'dimensions, or as a comma-separated list of the size for each '
+        'voxel dimension.',
+        xor=['image_size', 'scale_factor'],
+    )
+    scale_factor = traits.Tuple(
+        (traits.Float, traits.Float, traits.Float),
+        argstr='-scale %d,%d,%d',
+        mandatory=True,
+        desc='scale the image resolution by the supplied factor. This can be '
+        'specified either as a single value to be used for all '
+        'dimensions, or as a comma-separated list of scale factors for '
+        'each dimension.',
+        xor=['image_size', 'voxel_size'],
+    )
+    interp = traits.Enum(
+        'cubic',
+        'nearest',
+        'linear',
+        'sinc',
+        usedefault=True,
+        argstr='-interp %s',
+        desc='set the interpolation method to use when resizing (choices: '
+        'nearest, linear, cubic, sinc. Default: cubic).',
+    )
+    out_file = File(
+        argstr='%s',
+        name_template='%s_resized',
+        name_source=['in_file'],
+        keep_extension=True,
+        position=-1,
+        desc='the output resized DWI image',
+    )
+
+
+class MRResizeOutputSpec(TraitedSpec):
+    out_file = File(desc='the output resized DWI image', exists=True)
+
+
+class MRResize(MRTrix3Base):
+    """
+    Resize an image by defining the new image resolution, voxel size or a
+    scale factor. If the image is 4D, then only the first 3 dimensions can be
+    resized. Also, if the image is down-sampled, the appropriate smoothing is
+    automatically applied using Gaussian smoothing.
+    For more information, see
+    <https://mrtrix.readthedocs.io/en/latest/reference/commands/mrresize.html>
+    Example
+    -------
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> resize = mrt.MRResize()
+    >>> resize.inputs.in_file = 'dwi.mif'
+    >>> resize.inputs.voxel_size = (1, 1, 1)
+    >>> resize.cmdline                               # doctest: +ELLIPSIS
+    'mrresize -interp cubic -voxel 1,1,1 dwi.mif dwi_resized.mif'
+    >>> resize.run()                                 # doctest: +SKIP
+    """
+
+    _cmd = 'mrresize'
+    input_spec = MRResizeInputSpec
+    output_spec = MRResizeOutputSpec
diff --git a/dmriprep/utils/bids.py b/dmriprep/utils/bids.py
index 8a3f4269..ec02a1d2 100644
--- a/dmriprep/utils/bids.py
+++ b/dmriprep/utils/bids.py
@@ -8,7 +8,7 @@
 from bids import BIDSLayout
 
 
-def collect_data(bids_dir, participant_label, task=None, echo=None,
+def collect_data(bids_dir, participant_label, acq_id=None, concat_dwis=None,
                  bids_validate=True):
     """Replacement for niworkflows' version."""
     if isinstance(bids_dir, BIDSLayout):
@@ -19,7 +19,6 @@ def collect_data(bids_dir, participant_label, task=None, echo=None,
     queries = {
         'fmap': {'datatype': 'fmap'},
         'dwi': {'datatype': 'dwi', 'suffix': 'dwi'},
-        'bold': {'datatype': 'func', 'suffix': 'bold'},
         'sbref': {'datatype': 'func', 'suffix': 'sbref'},
         'flair': {'datatype': 'anat', 'suffix': 'FLAIR'},
         't2w': {'datatype': 'anat', 'suffix': 'T2w'},
@@ -27,17 +26,13 @@ def collect_data(bids_dir, participant_label, task=None, echo=None,
         'roi': {'datatype': 'anat', 'suffix': 'roi'},
     }
 
-    if task:
-        queries['bold']['task'] = task
-
-    if echo:
-        queries['bold']['echo'] = echo
-
     subj_data = {
         dtype: sorted(layout.get(return_type='file', subject=participant_label,
                                  extension=['nii', 'nii.gz'], **query))
         for dtype, query in queries.items()}
 
+    # subj_data['dwi'] = group_dwi(subj_data['dwi'], acq_id, concat_dwis)
+
     return subj_data, layout
 
 
@@ -60,10 +55,10 @@ def write_derivative_description(bids_dir, deriv_dir):
     }
 
     # Keys that can only be set by environment
-    if 'FMRIPREP_DOCKER_TAG' in os.environ:
-        desc['DockerHubContainerTag'] = os.environ['FMRIPREP_DOCKER_TAG']
-    if 'FMRIPREP_SINGULARITY_URL' in os.environ:
-        singularity_url = os.environ['FMRIPREP_SINGULARITY_URL']
+    if 'DMRIPREP_DOCKER_TAG' in os.environ:
+        desc['DockerHubContainerTag'] = os.environ['DMRIPREP_DOCKER_TAG']
+    if 'DMRIPREP_SINGULARITY_URL' in os.environ:
+        singularity_url = os.environ['DMRIPREP_SINGULARITY_URL']
         desc['SingularityContainerURL'] = singularity_url
 
         singularity_md5 = _get_shub_version(singularity_url)
@@ -88,7 +83,7 @@ def write_derivative_description(bids_dir, deriv_dir):
 
 
 def validate_input_dir(exec_env, bids_dir, participant_label):
-    # Ignore issues and warnings that should not influence FMRIPREP
+    # Ignore issues and warnings that should not influence DMRIPREP
     import tempfile
     import subprocess
     validator_config_dict = {
@@ -98,25 +93,20 @@ def validate_input_dir(exec_env, bids_dir, participant_label):
             "TSV_EQUAL_ROWS",
             "TSV_EMPTY_CELL",
             "TSV_IMPROPER_NA",
-            "VOLUME_COUNT_MISMATCH",
-            "BVAL_MULTIPLE_ROWS",
-            "BVEC_NUMBER_ROWS",
-            "DWI_MISSING_BVAL",
             "INCONSISTENT_SUBJECTS",
             "INCONSISTENT_PARAMETERS",
-            "BVEC_ROW_LENGTH",
-            "B_FILE",
             "PARTICIPANT_ID_COLUMN",
             "PARTICIPANT_ID_MISMATCH",
             "TASK_NAME_MUST_DEFINE",
             "PHENOTYPE_SUBJECTS_MISSING",
             "STIMULUS_FILE_MISSING",
-            "DWI_MISSING_BVEC",
+            "BOLD_NOT_4D",
             "EVENTS_TSV_MISSING",
             "TSV_IMPROPER_NA",
             "ACQTIME_FMT",
             "Participants age 89 or higher",
             "DATASET_DESCRIPTION_JSON_MISSING",
+            "TASK_NAME_CONTAIN_ILLEGAL_CHARACTER",
             "FILENAME_COLUMN",
             "WRONG_NEW_LINE",
             "MISSING_TSV_COLUMN_CHANNELS",
@@ -131,8 +121,6 @@ def validate_input_dir(exec_env, bids_dir, participant_label):
             "ACQTIME_FMT",
             "SUSPICIOUSLY_LONG_EVENT_DESIGN",
             "SUSPICIOUSLY_SHORT_EVENT_DESIGN",
-            "MALFORMED_BVEC",
-            "MALFORMED_BVAL",
             "MISSING_TSV_COLUMN_EEG_ELECTRODES",
             "MISSING_SESSION"
         ],
@@ -174,5 +162,12 @@ def validate_input_dir(exec_env, bids_dir, participant_label):
             print("bids-validator does not appear to be installed", file=sys.stderr)
 
 
+# def group_dwi(dwi_files, acq_id=None, concat_dwis=None):
+#
+#     all_dwis = []
+#
+#     return all_dwis
+
+
 def _get_shub_version(singularity_url):
     return NotImplemented
diff --git a/dmriprep/utils/hemisphere.py b/dmriprep/utils/hemisphere.py
new file mode 100644
index 00000000..648f2a2f
--- /dev/null
+++ b/dmriprep/utils/hemisphere.py
@@ -0,0 +1,62 @@
+import itertools
+import logging
+
+import numpy as np
+
+
+def is_hemispherical(vecs):
+    """
+    Test whether all points on a unit sphere lie in the same hemisphere.
+
+    **Inputs**
+
+    vecs : numpy.ndarray
+        2D numpy array with shape (N, 3) where N is the number of points.
+        All points must lie on the unit sphere.
+
+    **Outputs**
+
+    is_hemi : bool
+        If True, one can find a hemisphere that contains all the points.
+        If False, then the points do not lie in any hemisphere
+    pole : numpy.ndarray
+        If `is_hemi == True`, then pole is the "central" pole of the
+        input vectors. Otherwise, pole is the zero vector.
+
+    **References**
+
+    https://rstudio-pubs-static.s3.amazonaws.com/27121_a22e51b47c544980bad594d5e0bb2d04.html  # noqa
+    """
+    if vecs.shape[1] != 3:
+        raise ValueError("Input vectors must be 3D vectors")
+    if not np.allclose(1, np.linalg.norm(vecs, axis=1)):
+        raise ValueError("Input vectors must be unit vectors")
+
+    # Generate all pairwise cross products
+    v0, v1 = zip(*[p for p in itertools.permutations(vecs, 2)])
+    cross_prods = np.cross(v0, v1)
+
+    # Normalize them
+    cross_prods /= np.linalg.norm(cross_prods, axis=1)[:, np.newaxis]
+
+    # `cross_prods` now contains all candidate vertex points for "the polygon"
+    # in the reference. "The polygon" is a subset. Find which points belong to
+    # the polygon using a dot product test with each of the original vectors
+    angles = np.arccos(np.dot(cross_prods, vecs.transpose()))
+
+    # And test whether it is orthogonal or less
+    dot_prod_test = angles <= np.pi / 2.0
+
+    # If there is at least one point that is orthogonal or less to each
+    # input vector, then the points lie on some hemisphere
+    is_hemi = len(vecs) in np.sum(dot_prod_test.astype(int), axis=1)
+
+    if is_hemi:
+        vertices = cross_prods[
+            np.sum(dot_prod_test.astype(int), axis=1) == len(vecs)
+        ]
+        pole = np.mean(vertices, axis=0)
+        pole /= np.linalg.norm(pole)
+    else:
+        pole = np.array([0.0, 0.0, 0.0])
+    return is_hemi, pole
diff --git a/dmriprep/utils/sentry.py b/dmriprep/utils/sentry.py
new file mode 100644
index 00000000..c2895fa6
--- /dev/null
+++ b/dmriprep/utils/sentry.py
@@ -0,0 +1,208 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""Stripped out routines for Sentry"""
+import os
+from pathlib import Path
+import re
+from niworkflows.utils.misc import read_crashfile
+import sentry_sdk
+
+CHUNK_SIZE = 16384
+# Group common events with pre specified fingerprints
+KNOWN_ERRORS = {
+    'permission-denied': [
+        "PermissionError: [Errno 13] Permission denied"
+    ],
+    'memory-error': [
+        "MemoryError",
+        "Cannot allocate memory",
+        "Return code: 134",
+    ],
+    'reconall-already-running': [
+        "ERROR: it appears that recon-all is already running"
+    ],
+    'no-disk-space': [
+        "[Errno 28] No space left on device",
+        "[Errno 122] Disk quota exceeded"
+    ],
+    'segfault': [
+        "Segmentation Fault",
+        "Segfault",
+        "Return code: 139",
+    ],
+    'potential-race-condition': [
+        "[Errno 39] Directory not empty",
+        "_unfinished.json",
+    ],
+    'keyboard-interrupt': [
+        "KeyboardInterrupt",
+    ],
+}
+
+
+def start_ping(run_uuid, npart):
+    with sentry_sdk.configure_scope() as scope:
+        if run_uuid:
+            scope.set_tag('run_uuid', run_uuid)
+        scope.set_tag('npart', npart)
+    sentry_sdk.add_breadcrumb(message='dMRIPrep started', level='info')
+    sentry_sdk.capture_message('dMRIPrep started', level='info')
+
+
+def sentry_setup(opts, exec_env):
+    from os import cpu_count
+    import psutil
+    import hashlib
+    from ..__about__ import __version__
+
+    environment = "prod"
+    release = __version__
+    if not __version__:
+        environment = "dev"
+        release = "dev"
+    elif int(os.getenv('DMRIPREP_DEV', '0')) or ('+' in __version__):
+        environment = "dev"
+
+    sentry_sdk.init("https://d5a16b0c38d84d1584dfc93b9fb1ade6@sentry.io/1137693",
+                    release=release,
+                    environment=environment,
+                    before_send=before_send)
+    with sentry_sdk.configure_scope() as scope:
+        scope.set_tag('exec_env', exec_env)
+
+        if exec_env == 'dmriprep-docker':
+            scope.set_tag('docker_version', os.getenv('DOCKER_VERSION_8395080871'))
+
+        dset_desc_path = opts.bids_dir / 'dataset_description.json'
+        if dset_desc_path.exists():
+            desc_content = dset_desc_path.read_bytes()
+            scope.set_tag('dset_desc_sha256', hashlib.sha256(desc_content).hexdigest())
+
+        free_mem_at_start = round(psutil.virtual_memory().free / 1024**3, 1)
+        scope.set_tag('free_mem_at_start', free_mem_at_start)
+        scope.set_tag('cpu_count', cpu_count())
+
+        # Memory policy may have a large effect on types of errors experienced
+        overcommit_memory = Path('/proc/sys/vm/overcommit_memory')
+        if overcommit_memory.exists():
+            policy = {'0': 'heuristic',
+                      '1': 'always',
+                      '2': 'never'}.get(overcommit_memory.read_text().strip(), 'unknown')
+            scope.set_tag('overcommit_memory', policy)
+            if policy == 'never':
+                overcommit_kbytes = Path('/proc/sys/vm/overcommit_memory')
+                kb = overcommit_kbytes.read_text().strip()
+                if kb != '0':
+                    limit = '{}kB'.format(kb)
+                else:
+                    overcommit_ratio = Path('/proc/sys/vm/overcommit_ratio')
+                    limit = '{}%'.format(overcommit_ratio.read_text().strip())
+                scope.set_tag('overcommit_limit', limit)
+            else:
+                scope.set_tag('overcommit_limit', 'n/a')
+        else:
+            scope.set_tag('overcommit_memory', 'n/a')
+            scope.set_tag('overcommit_limit', 'n/a')
+
+        for k, v in vars(opts).items():
+            scope.set_tag(k, v)
+
+
+def process_crashfile(crashfile):
+    """Parse the contents of a crashfile and submit sentry messages"""
+    crash_info = read_crashfile(str(crashfile))
+    with sentry_sdk.push_scope() as scope:
+        scope.level = 'fatal'
+
+        # Extract node name
+        node_name = crash_info.pop('node').split('.')[-1]
+        scope.set_tag("node_name", node_name)
+
+        # Massage the traceback, extract the gist
+        traceback = crash_info.pop('traceback')
+        # last line is probably most informative summary
+        gist = traceback.splitlines()[-1]
+        exception_text_start = 1
+        for line in traceback.splitlines()[1:]:
+            if not line[0].isspace():
+                break
+            exception_text_start += 1
+
+        exception_text = '\n'.join(
+            traceback.splitlines()[exception_text_start:])
+
+        # Extract inputs, if present
+        inputs = crash_info.pop('inputs', None)
+        if inputs:
+            scope.set_extra('inputs', dict(inputs))
+
+        # Extract any other possible metadata in the crash file
+        for k, v in crash_info.items():
+            strv = list(_chunks(str(v)))
+            if len(strv) == 1:
+                scope.set_extra(k, strv[0])
+            else:
+                for i, chunk in enumerate(strv):
+                    scope.set_extra('%s_%02d' % (k, i), chunk)
+
+        fingerprint = ''
+        issue_title = '{}: {}'.format(node_name, gist)
+        for new_fingerprint, error_snippets in KNOWN_ERRORS.items():
+            for error_snippet in error_snippets:
+                if error_snippet in traceback:
+                    fingerprint = new_fingerprint
+                    issue_title = new_fingerprint
+                    break
+            if fingerprint:
+                break
+
+        message = issue_title + '\n\n'
+        message += exception_text[-(8192 - len(message)):]
+        if fingerprint:
+            sentry_sdk.add_breadcrumb(message=fingerprint, level='fatal')
+        else:
+            # remove file paths
+            fingerprint = re.sub(r"(/[^/ ]*)+/?", '', message)
+            # remove words containing numbers
+            fingerprint = re.sub(r"([a-zA-Z]*[0-9]+[a-zA-Z]*)+", '', fingerprint)
+            # adding the return code if it exists
+            for line in message.splitlines():
+                if line.startswith("Return code"):
+                    fingerprint += line
+                    break
+
+        scope.fingerprint = [fingerprint]
+        sentry_sdk.capture_message(message, 'fatal')
+
+
+def before_send(event, hints):
+    # Filtering log messages about crashed nodes
+    if 'logentry' in event and 'message' in event['logentry']:
+        msg = event['logentry']['message']
+        if msg.startswith("could not run node:"):
+            return None
+        if msg.startswith("Saving crash info to "):
+            return None
+        if re.match("Node .+ failed to run on host .+", msg):
+            return None
+
+    if 'breadcrumbs' in event and isinstance(event['breadcrumbs'], list):
+        fingerprints_to_propagate = ['no-disk-space', 'memory-error', 'permission-denied',
+                                     'keyboard-interrupt']
+        for bc in event['breadcrumbs']:
+            msg = bc.get('message', 'empty-msg')
+            if msg in fingerprints_to_propagate:
+                event['fingerprint'] = [msg]
+                break
+
+    return event
+
+
+def _chunks(string, length=CHUNK_SIZE):
+    """
+    Splits a string into smaller chunks
+    >>> list(_chunks('some longer string.', length=3))
+    ['som', 'e l', 'ong', 'er ', 'str', 'ing', '.']
+    """
+    return (string[i:i + length]
+            for i in range(0, len(string), length))
diff --git a/dmriprep/workflows/anatomical.py b/dmriprep/workflows/anatomical.py
new file mode 100644
index 00000000..7101b8bf
--- /dev/null
+++ b/dmriprep/workflows/anatomical.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+
+"""
+Orchestrating the T1w preprocessing workflows
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_anat_preproc_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+
+from niworkflows.anat.ants import init_brain_extraction_wf
+
+
+def init_anat_preproc_wf():
+    """
+    """
+
+    wf = pe.Workflow(name='init_anat_preproc_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=['T1w_file']), name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']), name='outputnode')
+
+    t1_skullstrip = init_brain_extraction_wf()
+
+    to_list = lambda x: [x]
+
+    wf.connect([
+        (inputnode, t1_skullstrip, [(('t1_file', to_list), 'inputnode.in_files')]),
+        (t1_skullstrip, outputnode, [('outputnode.out_file', 'out_file')])
+    ])
+
+    return wf
diff --git a/dmriprep/workflows/base.py b/dmriprep/workflows/base.py
old mode 100755
new mode 100644
index 2474204f..d29e596b
--- a/dmriprep/workflows/base.py
+++ b/dmriprep/workflows/base.py
@@ -29,12 +29,18 @@
 from ..interfaces.reports import SubjectSummary, AboutSummary
 from ..utils.bids import collect_data
 from ..__about__ import __version__
-# from .dwi import init_dwi_preproc_wf
+from .dwi import init_dwi_preproc_wf
 
 
 def init_dmriprep_wf(
+    acq_id,
+    acqp_file,
     anat_only,
+    b0_thresh,
+    concat_dwis,
     debug,
+    eddy_config,
+    fmap_bspline,
     force_syn,
     freesurfer,
     hires,
@@ -44,6 +50,7 @@ def init_dmriprep_wf(
     low_mem,
     omp_nthreads,
     output_dir,
+    output_resolution,
     output_spaces,
     run_uuid,
     skull_strip_fixed_seed,
@@ -69,8 +76,14 @@ def init_dmriprep_wf(
         from dmriprep.workflows.base import init_dmriprep_wf
         os.environ['FREESURFER_HOME'] = os.getcwd()
         wf = init_dmriprep_wf(
+            acq_id=[],
+            acqp_file='',
             anat_only=False,
+            b0_thresh=5,
+            concat_dwis=[],
             debug=False,
+            eddy_config='',
+            fmap_bspline=False,
             force_syn=True,
             freesurfer=True,
             hires=True,
@@ -80,6 +93,7 @@ def init_dmriprep_wf(
             low_mem=False,
             omp_nthreads=1,
             output_dir='.',
+            output_resolution=1,
             output_spaces=OrderedDict([
                 ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}),
                 ('T1w', {}), ('fsnative', {})]),
@@ -94,10 +108,22 @@ def init_dmriprep_wf(
 
     Parameters
     ----------
+        acq_id : list
+            List of acquisition labels of dwi scans to preprocess
+        acqp_file : str
+            Custom acquisition parameter file for topup/eddy
         anat_only : bool
             Disable diffusion MRI workflows
+        b0_thresh : int
+            Threshold of bvals to be classified as b0s
+        concat_dwis : list
+            List of acquisition labels of dwi scans to concatenate
         debug : bool
             Enable debugging outputs
+        eddy_config : str
+            Path to json file of eddy input arguments
+        fmap_bspline : bool
+            Fit a B-Spline using least-squares (experimental)
         force_syn : bool
             **Temporary**: Always run SyN-based SDC
         freesurfer : bool
@@ -105,7 +131,7 @@ def init_dmriprep_wf(
         hires : bool
             Enable sub-millimeter preprocessing in FreeSurfer
         ignore : list
-            Preprocessing steps to skip (may include "slicetiming", "fieldmaps")
+            Preprocessing steps to skip (may include "denoising", "unringing", "fieldmaps")
         layout : BIDSLayout object
             BIDS dataset layout
         longitudinal : bool
@@ -117,6 +143,8 @@ def init_dmriprep_wf(
             Maximum number of threads an individual process may use
         output_dir : str
             Directory in which to save derivatives
+        output_resolution : float
+            Output resolution of dwi scan
         output_spaces : OrderedDict
             Ordered dictionary where keys are TemplateFlow ID strings (e.g., ``MNI152Lin``,
             ``MNI152NLin6Asym``, ``MNI152NLin2009cAsym``, or ``fsLR``) strings designating
@@ -157,8 +185,14 @@ def init_dmriprep_wf(
     reportlets_dir = os.path.join(work_dir, 'reportlets')
     for subject_id in subject_list:
         single_subject_wf = init_single_subject_wf(
+            acq_id=acq_id,
+            acqp_file=acqp_file,
             anat_only=anat_only,
+            b0_thresh=b0_thresh,
+            concat_dwis=concat_dwis,
             debug=debug,
+            eddy_config=eddy_config,
+            fmap_bspline=fmap_bspline,
             force_syn=force_syn,
             freesurfer=freesurfer,
             hires=hires,
@@ -169,6 +203,7 @@ def init_dmriprep_wf(
             name="single_subject_" + subject_id + "_wf",
             omp_nthreads=omp_nthreads,
             output_dir=output_dir,
+            output_resolution=output_resolution,
             output_spaces=output_spaces,
             reportlets_dir=reportlets_dir,
             skull_strip_fixed_seed=skull_strip_fixed_seed,
@@ -192,8 +227,14 @@ def init_dmriprep_wf(
 
 
 def init_single_subject_wf(
+    acq_id,
+    acqp_file,
     anat_only,
+    b0_thresh,
+    concat_dwis,
     debug,
+    eddy_config,
+    fmap_bspline,
     force_syn,
     freesurfer,
     hires,
@@ -204,6 +245,7 @@ def init_single_subject_wf(
     name,
     omp_nthreads,
     output_dir,
+    output_resolution,
     output_spaces,
     reportlets_dir,
     skull_strip_fixed_seed,
@@ -230,8 +272,14 @@ def init_single_subject_wf(
         from collections import namedtuple, OrderedDict
         BIDSLayout = namedtuple('BIDSLayout', ['root'])
         wf = init_single_subject_wf(
+            acq_id=[],
+            acqp_file='',
             anat_only=False,
+            b0_thresh=5,
+            concat_dwis=[],
             debug=False,
+            eddy_config='',
+            fmap_bspline=False,
             force_syn=True,
             freesurfer=True,
             hires=True,
@@ -242,6 +290,7 @@ def init_single_subject_wf(
             name='single_subject_wf',
             omp_nthreads=1,
             output_dir='.',
+            output_resolution=1,
             output_spaces=OrderedDict([
                 ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}),
                 ('T1w', {}), ('fsnative', {})]),
@@ -255,10 +304,22 @@ def init_single_subject_wf(
 
     Parameters
     ----------
+        acq_id : list
+            List of acquisition labels of dwi scans to preprocess
+        acqp_file : str
+            Custom acquisition parameter file for topup/eddy
         anat_only : bool
             Disable diffusion MRI workflows
+        b0_thresh : int
+            Threshold of bvals to be classified as b0s
+        concat_dwis : list
+            List of acquisition labels of dwi scans to concatenate
         debug : bool
             Enable debugging outputs
+        eddy_config : str
+            Path to json file of eddy input arguments
+        fmap_bspline : bool
+            Fit a B-Spline using least-squares (experimental)
         force_syn : bool
             **Temporary**: Always run SyN-based SDC
         freesurfer : bool
@@ -280,6 +341,8 @@ def init_single_subject_wf(
             Maximum number of threads an individual process may use
         output_dir : str
             Directory in which to save derivatives
+        output_resolution : float
+            Output resolution of dwi scan
         output_spaces : OrderedDict
             Ordered dictionary where keys are TemplateFlow ID strings (e.g., ``MNI152Lin``,
             ``MNI152NLin6Asym``, ``MNI152NLin2009cAsym``, or ``fsLR``) strings designating
@@ -314,7 +377,7 @@ def init_single_subject_wf(
         # for documentation purposes
         subject_data = {
             't1w': ['/completely/made/up/path/sub-01_T1w.nii.gz'],
-            'dwi': ['/completely/made/up/path/sub-01_task-nback_bold.nii.gz']
+            'dwi': ['/completely/made/up/path/sub-01_dwi.nii.gz']
         }
     else:
         subject_data = collect_data(layout, subject_id)[0]
@@ -437,59 +500,52 @@ def init_single_subject_wf(
     if anat_only:
         return workflow
 
-    # for dwi_file in subject_data['dwi']:
-    #     dwi_preproc_wf = init_dwi_preproc_wf(
-    #         aroma_melodic_dim=aroma_melodic_dim,
-    #         bold2t1w_dof=bold2t1w_dof,
-    #         bold_file=bold_file,
-    #         cifti_output=cifti_output,
-    #         debug=debug,
-    #         dummy_scans=dummy_scans,
-    #         err_on_aroma_warn=err_on_aroma_warn,
-    #         fmap_bspline=fmap_bspline,
-    #         fmap_demean=fmap_demean,
-    #         force_syn=force_syn,
-    #         freesurfer=freesurfer,
-    #         ignore=ignore,
-    #         layout=layout,
-    #         low_mem=low_mem,
-    #         medial_surface_nan=medial_surface_nan,
-    #         num_bold=len(subject_data['bold']),
-    #         omp_nthreads=omp_nthreads,
-    #         output_dir=output_dir,
-    #         output_spaces=output_spaces,
-    #         reportlets_dir=reportlets_dir,
-    #         regressors_all_comps=regressors_all_comps,
-    #         regressors_fd_th=regressors_fd_th,
-    #         regressors_dvars_th=regressors_dvars_th,
-    #         t2s_coreg=t2s_coreg,
-    #         use_aroma=use_aroma,
-    #         use_syn=use_syn,
-    #     )
-
-    #     workflow.connect([
-    #         (anat_preproc_wf, dwi_preproc_wf,
-    #          [(('outputnode.t1_preproc', _pop), 'inputnode.t1_preproc'),
-    #           ('outputnode.t1_brain', 'inputnode.t1_brain'),
-    #           ('outputnode.t1_mask', 'inputnode.t1_mask'),
-    #           ('outputnode.t1_seg', 'inputnode.t1_seg'),
-    #           ('outputnode.t1_aseg', 'inputnode.t1_aseg'),
-    #           ('outputnode.t1_aparc', 'inputnode.t1_aparc'),
-    #           ('outputnode.t1_tpms', 'inputnode.t1_tpms'),
-    #           ('outputnode.template', 'inputnode.template'),
-    #           ('outputnode.forward_transform', 'inputnode.anat2std_xfm'),
-    #           ('outputnode.reverse_transform', 'inputnode.std2anat_xfm'),
-    #           ('outputnode.joint_template', 'inputnode.joint_template'),
-    #           ('outputnode.joint_forward_transform', 'inputnode.joint_anat2std_xfm'),
-    #           ('outputnode.joint_reverse_transform', 'inputnode.joint_std2anat_xfm'),
-    #           # Undefined if --no-freesurfer, but this is safe
-    #           ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
-    #           ('outputnode.subject_id', 'inputnode.subject_id'),
-    #           ('outputnode.t1_2_fsnative_forward_transform',
-    #            'inputnode.t1_2_fsnative_forward_transform'),
-    #           ('outputnode.t1_2_fsnative_reverse_transform',
-    #            'inputnode.t1_2_fsnative_reverse_transform')]),
-    #     ])
+    for dwi_file in subject_data['dwi']:
+        dwi_preproc_wf = init_dwi_preproc_wf(
+            acqp_file=acqp_file,
+            b0_thresh=b0_thresh,
+            debug=debug,
+            dwi_file=dwi_file,
+            eddy_config=eddy_config,
+            fmap_bspline=fmap_bspline,
+            # fmap_demean=fmap_demean,
+            force_syn=force_syn,
+            freesurfer=freesurfer,
+            ignore=ignore,
+            layout=layout,
+            low_mem=low_mem,
+            num_dwi=len(subject_data['dwi']),
+            omp_nthreads=omp_nthreads,
+            output_dir=output_dir,
+            output_resolution=output_resolution,
+            output_spaces=output_spaces,
+            reportlets_dir=reportlets_dir,
+            use_syn=use_syn,
+        )
+
+        workflow.connect([
+            (anat_preproc_wf, dwi_preproc_wf,
+             [#(('outputnode.t1_preproc', _pop), 'inputnode.t1_preproc'),
+              ('outputnode.t1_brain', 'inputnode.t1_brain'),
+              # ('outputnode.t1_mask', 'inputnode.t1_mask'),
+              # ('outputnode.t1_seg', 'inputnode.t1_seg'),
+              # ('outputnode.t1_aseg', 'inputnode.t1_aseg'),
+              # ('outputnode.t1_aparc', 'inputnode.t1_aparc'),
+              # ('outputnode.t1_tpms', 'inputnode.t1_tpms'),
+              ('outputnode.template', 'inputnode.template'),
+              # ('outputnode.forward_transform', 'inputnode.anat2std_xfm'),
+              ('outputnode.reverse_transform', 'inputnode.std2anat_xfm'),
+              # ('outputnode.joint_template', 'inputnode.joint_template'),
+              # ('outputnode.joint_forward_transform', 'inputnode.joint_anat2std_xfm'),
+              # ('outputnode.joint_reverse_transform', 'inputnode.joint_std2anat_xfm'),
+              # Undefined if --no-freesurfer, but this is safe
+              ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
+              ('outputnode.subject_id', 'inputnode.subject_id')]),
+              # ('outputnode.t1_2_fsnative_forward_transform',
+              #  'inputnode.t1_2_fsnative_forward_transform'),
+              # ('outputnode.t1_2_fsnative_reverse_transform',
+              #  'inputnode.t1_2_fsnative_reverse_transform')]),
+        ])
 
     return workflow
 
diff --git a/dmriprep/workflows/dwi/__init__.py b/dmriprep/workflows/dwi/__init__.py
index e69de29b..401e8771 100644
--- a/dmriprep/workflows/dwi/__init__.py
+++ b/dmriprep/workflows/dwi/__init__.py
@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+
+"""
+
+Pre-processing dMRI workflows
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: dmriprep.workflows.dwi.base
+.. automodule:: dmriprep.workflows.dwi.util
+.. automodule:: dmriprep.workflows.dwi.artifacts
+.. automodule:: dmriprep.workflows.dwi.eddy
+.. automodule:: dmriprep.workflows.dwi.tensor
+.. automodule:: dmriprep.workflows.dwi.outputs
+
+"""
+
+from .base import init_dwi_preproc_wf
+from .util import init_dwi_concat_wf
+from .artifacts import init_dwi_artifacts_wf
+from .eddy import init_dwi_eddy_wf
+from .tensor import init_dwi_tensor_wf
+from .outputs import init_dwi_derivatives_wf
+
+__all__ = [
+    'init_dwi_preproc_wf',
+    'init_dwi_concat_wf',
+    'init_dwi_artifacts_wf',
+    'init_dwi_eddy_wf',
+    'init_dwi_tensor_wf',
+    'init_dwi_derivatives_wf'
+]
diff --git a/dmriprep/workflows/dwi/artifacts.py b/dmriprep/workflows/dwi/artifacts.py
new file mode 100644
index 00000000..c7917244
--- /dev/null
+++ b/dmriprep/workflows/dwi/artifacts.py
@@ -0,0 +1,132 @@
+# -*- coding: utf-8 -*-
+
+"""
+Artefact removal and resizing
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_dwi_artifacts_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+
+from dmriprep.due import due, BibTeX
+from dmriprep.interfaces import mrtrix3
+
+
+@due.dcite(BibTeX(
+    """
+    @article{denoise1,
+    title = {Denoising of diffusion MRI using random matrix theory},
+    volume = {142},
+    doi = {10.1016/j.neuroimage.2016.08.016},
+    journal = {NeuroImage},
+    author = {Veraart, J.; Novikov, D.S.; Christiaens, D.; Ades-aron, B.; Sijbers, J. & Fieremans, E.},
+    year = {2016},
+    pages = {394-406}}
+    """),
+    description='Denoising algorithm')
+@due.dcite(BibTeX(
+    """
+    @article{denoise2,
+    title = {Diffusion MRI noise mapping using random matrix theory},
+    volume = {76},
+    doi = {10.1002/mrm.26059},
+    journal = {Magnetic Resonance in Medicine},
+    author = {Veraart, J.; Fieremans, E. & Novikov, D.S.},
+    year = {2016},
+    pages = {1582-1593}}
+    """),
+    description='Denoising algorithm')
+@due.dcite(BibTeX(
+    """
+    @article{unring,
+    title = {Gibbs-ringing artifact removal based on local subvoxel-shifts},
+    volume = {76},
+    doi = {10.1002/mrm.26054},
+    journal = {Magnetic Resonance in Medicine},
+    author = {Kellner, E; Dhital, B; Kiselev, V.G & Reisert, M.},
+    year = {2016},
+    pages = {1574–1581}}
+    """),
+    description='Unringing algorithm')
+def init_dwi_artifacts_wf(ignore, output_resolution):
+    """
+    This workflow performs denoising and unringing on the input dwi image.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from dmriprep.workflows.dwi import init_dwi_artifacts_wf
+        wf = init_dwi_artifacts_wf(ignore=[], output_resolution=(1, 1, 1))
+
+    **Parameters**
+
+        ignore : :obj:`list`
+            List of artefact removal steps to skip (default: [])
+
+    **Inputs**
+
+        dwi_file : :obj:`str`
+            dwi NIfTI file
+
+    **Outputs**
+
+        out_file : :obj:`str`
+            dwi NIfTI file after artefact removal
+
+    """
+
+    wf = pe.Workflow(name='dwi_artifacts_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=['dwi_file']), name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']), name='outputnode')
+
+    dwibuffer = pe.Node(niu.IdentityInterface(fields=['dwi_file']), name='dwibuffer')
+
+    denoise = pe.Node(mrtrix3.DWIDenoise(), name='denoise')
+
+    unring = pe.Node(mrtrix3.MRDeGibbs(), name='unring')
+
+    resize = pe.Node(mrtrix3.MRResize(), name='resize')
+
+    if ignore == ['denoising']:
+        wf.connect([
+            (inputnode, unring, [('dwi_file', 'in_file')]),
+            (unring, dwibuffer, [('out_file', 'dwi_file')])
+        ])
+
+    elif ignore == ['unringing']:
+        wf.connect([
+            (inputnode, denoise, [('dwi_file', 'in_file')]),
+            (denoise, dwibuffer, [('out_file', 'dwi_file')])
+        ])
+
+    elif ['denoising', 'unringing'] in ignore:
+        wf.connect([
+            (inputnode, dwibuffer, 'dwi_file', 'dwi_file')
+        ])
+
+    else:
+        wf.connect([
+            (inputnode, denoise, [('dwi_file', 'in_file')]),
+            (denoise, unring, [('out_file', 'in_file')]),
+            (unring, dwibuffer, [('out_file', 'dwi_file')])
+        ])
+
+    if output_resolution:
+        resize.inputs.voxel_size = output_resolution
+        wf.connect([
+            (dwibuffer, resize, [('dwi_file', 'in_file')]),
+            (resize, outputnode, [('out_file', 'out_file')])
+        ])
+
+    else:
+        wf.connect([
+            (dwibuffer, outputnode, [('dwi_file', 'out_file')])
+        ])
+
+    return wf
diff --git a/dmriprep/workflows/dwi/base.py b/dmriprep/workflows/dwi/base.py
new file mode 100755
index 00000000..af12eab3
--- /dev/null
+++ b/dmriprep/workflows/dwi/base.py
@@ -0,0 +1,437 @@
+# -*- coding: utf-8 -*-
+
+"""
+Orchestrating the dwi preprocessing workflows
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_dwi_preproc_wf
+
+"""
+
+from nipype import logging
+from nipype.pipeline import engine as pe
+from nipype.interfaces import fsl, mrtrix3, utility as niu
+
+from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+from .util import init_dwi_concat_wf
+from .artifacts import init_dwi_artifacts_wf
+from .eddy import init_dwi_eddy_wf
+from .tensor import init_dwi_tensor_wf
+from ..fieldmap.base import init_sdc_wf
+
+FMAP_PRIORITY = {'epi': 0, 'fieldmap': 1, 'phasediff': 2, 'phase': 3, 'syn': 4}
+
+DEFAULT_MEMORY_MIN_GB = 0.01
+LOGGER = logging.getLogger('nipype.workflow')
+
+
+def init_dwi_preproc_wf(
+    acqp_file,
+    b0_thresh,
+    debug,
+    dwi_file,
+    eddy_config,
+    fmap_bspline,
+    force_syn,
+    freesurfer,
+    ignore,
+    layout,
+    low_mem,
+    num_dwi,
+    omp_nthreads,
+    output_dir,
+    output_resolution,
+    output_spaces,
+    reportlets_dir,
+    use_syn,
+):
+    """
+    This workflow controls the diffusion preprocessing stages of dMRIPrep.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from collections import namedtuple
+        from dmriprep.workflows.dwi import init_dwi_preproc_wf
+        BIDSLayout = namedtuple('BIDSLayout', ['root'])
+        wf = init_dwi_preproc_wf(
+            acqp_file='',
+            b0_thresh=5,
+            debug=False,
+            dwi_file='/madeup/path/sub-01_dwi.nii.gz',
+            eddy_config='',
+            fmap_bspline=False,
+            force_syn=False,
+            freesurfer=True,
+            ignore=[],
+            layout=BIDSLayout('.'),
+            low_mem=False,
+            num_dwi=1,
+            omp_nthreads=1,
+            output_dir='.',
+            output_resolution=(1, 1, 1),
+            output_spaces=OrderedDict([('dwi', {})]),
+            reportlets_dir='.',
+            use_syn=True
+        )
+
+    **Subworkflows**
+
+        * :py:func:`~dmriprep.workflows.dwi.util.init_dwi_concat_wf`
+        * :py:func:`~dmriprep.workflows.dwi.artifacts.init_dwi_artifacts_wf`
+        * :py:func:`~dmriprep.workflows.dwi.eddy.init_dwi_eddy_wf`
+        * :py:func:`~dmriprep.workflows.dwi.tensor.init_dwi_tensor_wf`
+
+    """
+
+    wf_name = _get_wf_name(dwi_file)
+
+    ref_file = dwi_file
+
+    if isinstance(dwi_file, list):
+        ref_file = dwi_file[0]
+
+    # For doc building purposes
+    if not hasattr(layout, 'parse_file_entities'):
+        metadata = {
+            'PhaseEncodingDirection': 'j-',
+            'TotalReadoutTime': 0.05
+        }
+        fmaps = [{
+            'suffix': 'phasediff',
+            'phasediff': 'sub-01/fmap/sub-01_phasediff.nii.gz',
+            'magnitude1': 'sub-01/fmap/sub-01_magnitude1.nii.gz',
+            'magnitude2': 'sub-01/fmap/sub-01_magnitude2.nii.gz'
+        }]
+    else:
+        metadata = layout.get_metadata(ref_file)
+        # Find fieldmaps. Options: (epi|fieldmap|phasediff|phase1|phase2|syn)
+        fmaps = []
+        if 'fieldmaps' not in ignore:
+            for fmap in layout.get_fieldmap(ref_file, return_list=True):
+                fmap['metadata'] = layout.get_metadata(fmap[fmap['suffix']])
+                fmaps.append(fmap)
+
+        if force_syn or (use_syn and not fmaps):
+            fmaps.append({'suffix': 'syn'})
+
+        if fmaps:
+            fmaps.sort(key=lambda fmap: FMAP_PRIORITY[fmap['suffix']])
+            fmap = fmaps[0]
+            if fmap['suffix'] == 'epi':
+                sdc_method = 'topup'
+            if any(s in fmap['suffix'] for s in ['fieldmap', 'phasediff', 'phase1']):
+                sdc_method = 'fieldmap'
+
+    dwi_wf = Workflow(name=wf_name)
+    dwi_wf.__desc__ = """
+
+Diffusion data preprocessing
+
+: For each of the {num_dwi} diffusion runs found per subject (across all
+sessions), the following preprocessing was performed.
+""".format(num_dwi=num_dwi)
+
+    inputnode = pe.Node(niu.IdentityInterface(
+        fields=['dwi_file', 'bvec_file', 'bval_file',
+                'subjects_dir', 'subject_id',
+                't1_brain', 'template', 'std2anat_xfm']),
+        name='inputnode')
+
+    if isinstance(dwi_file, list):
+        dwi_concat_wf = init_dwi_concat_wf(layout)
+        dwi_concat_wf.inputs.inputnode.ref_file = ref_file
+        dwi_concat_wf.inputs.inputnode.dwi_list = dwi_file
+
+        dwi_wf.connect([
+            (dwi_concat_wf, inputnode, [('outputnode.dwi_file' 'dwi_file'),
+                                        ('outputnode.bvec_file', 'bvec_file'),
+                                        ('outputnode.bval_file', 'bval_file')])
+        ])
+    else:
+        inputnode.inputs.dwi_file = dwi_file
+        inputnode.inputs.bvec_file = layout.get_bvec(dwi_file)
+        inputnode.inputs.bval_file = layout.get_bval(dwi_file)
+
+    outputnode = pe.Node(niu.IdentityInterface(
+        fields=['out_dwi', 'out_bval', 'out_bvec', 'index', 'acq_params',
+                'out_mask', 'out_b0_pre', 'out_b0_mask_pre',
+                'out_eddy_quad_json', 'out_eddy_quad_pdf',
+                'out_dtifit_FA', 'out_dtifit_V1', 'out_dtifit_sse',
+                'out_noise']),
+        name='outputnode')
+
+    def gen_index(in_file):
+        """
+        Create an index file that tells `eddy` which line/lines in the
+        `acqparams.txt` file are relevant for the dwi data passed into `eddy`
+
+        **Inputs**
+
+        in_file
+            The bval file.
+
+        **Outputs**
+
+        out_file
+            The output index file.
+
+        """
+
+        import os
+        import numpy as np
+        import nibabel as nib
+        from nipype.utils import NUMPY_MMAP
+        from nipype.utils.filemanip import fname_presuffix
+
+        out_file = fname_presuffix(
+            in_file, suffix='_index.txt', newpath=os.path.abspath('.'), use_ext=False
+        )
+        vols = nib.load(in_file, mmap=NUMPY_MMAP).get_data().shape[-1]
+        index_lines = np.ones((vols,))
+        index_lines_reshape = index_lines.reshape(1, index_lines.shape[0])
+        np.savetxt(out_file, index_lines_reshape, fmt='%i')
+        return out_file
+
+    gen_idx = pe.Node(
+        niu.Function(
+            input_names=['in_file'], output_names=['out_file'], function=gen_index),
+        name='gen_index')
+
+    def gen_acqparams(in_file, metadata):
+        """
+        Create an acquisition parameters file for ``eddy``
+
+        **Inputs**
+
+        in_file
+            The dwi file
+        metadata
+            The BIDS metadata of the dwi file
+
+        **Outputs**
+
+        out_file
+            The output acquisition parameters file
+
+        """
+
+        import os
+        from nipype.utils.filemanip import fname_presuffix
+
+        out_file = fname_presuffix(
+            in_file,
+            suffix='_acqparams.txt',
+            newpath=os.path.abspath('.'),
+            use_ext=False,
+        )
+
+        acq_param_dict = {
+            'j': '0 1 0 %.7f',
+            'j-': '0 -1 0 %.7f',
+            'i': '1 0 0 %.7f',
+            'i-': '-1 0 0 %.7f',
+            'k': '0 0 1 %.7f',
+            'k-': '0 0 -1 %.7f',
+        }
+
+        pe_dir = metadata.get('PhaseEncodingDirection')
+        total_readout = metadata.get('TotalReadoutTime')
+        acq_param_lines = acq_param_dict[pe_dir] % total_readout
+
+        with open(out_file, 'w') as f:
+            f.write(acq_param_lines)
+
+        return out_file
+
+    acqp = pe.Node(
+        niu.Function(
+            input_names=['in_file', 'metadata'],
+            output_names=['out_file'],
+            function=gen_acqparams),
+        name='acqp')
+    acqp.inputs.metadata = metadata
+
+    dwi_artifacts_wf = init_dwi_artifacts_wf(ignore, output_resolution)
+
+    def b0_average(in_dwi, in_bval, b0_thresh, out_file=None):
+        """
+        Averages the *b0* volumes from a DWI dataset.
+        As current dMRI data are being acquired with all b-values > 0.0,
+        the *lowb* volumes are selected by specifying the parameter b0_thresh.
+        .. warning:: *b0* should be already registered (head motion artifact
+        should be corrected).
+        """
+        import os
+        import numpy as np
+        import nibabel as nib
+        from nipype.utils import NUMPY_MMAP
+        from nipype.utils.filemanip import fname_presuffix
+
+        if out_file is None:
+            out_file = fname_presuffix(
+                in_dwi, suffix='_avg_b0', newpath=os.path.abspath('.')
+            )
+
+        imgs = np.array(nib.four_to_three(nib.load(in_dwi, mmap=NUMPY_MMAP)))
+        bval = np.loadtxt(in_bval)
+        index = np.argwhere(bval <= b0_thresh).flatten().tolist()
+
+        b0s = [im.get_data().astype(np.float32) for im in imgs[index]]
+        b0 = np.average(np.array(b0s), axis=0)
+
+        hdr = imgs[0].header.copy()
+        hdr.set_data_shape(b0.shape)
+        hdr.set_xyzt_units('mm')
+        hdr.set_data_dtype(np.float32)
+        nib.Nifti1Image(b0, imgs[0].affine, hdr).to_filename(out_file)
+        return out_file
+
+    avg_b0_0 = pe.Node(
+        niu.Function(
+            input_names=['in_dwi', 'in_bval', 'b0_thresh'],
+            output_names=['out_file'],
+            function=b0_average),
+        name='b0_avg_pre')
+
+    avg_b0_0.inputs.b0_thresh = b0_thresh
+
+    bet_dwi0 = pe.Node(
+        fsl.BET(frac=0.3, mask=True, robust=True), name='bet_dwi_pre'
+    )
+
+    dwi_sdc_wf = init_sdc_wf(
+        fmaps=fmaps,
+        metadata=metadata,
+        layout=layout,
+        ignore=ignore
+    )
+
+    dwi_eddy_wf = init_dwi_eddy_wf(omp_nthreads, sdc_method)
+
+    denoise_eddy = pe.Node(mrtrix3.DWIDenoise(), name='denoise_eddy')
+
+    eddy_quad = pe.Node(fsl.EddyQuad(verbose=True), name='eddy_quad')
+
+    get_path = lambda x: x.split('.nii.gz')[0].split('_fix')[0]
+    get_qc_path = lambda x: x.split('.nii.gz')[0] + '.qc'
+
+    fslroi = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name='fslroi')
+
+    bias_correct = pe.Node(mrtrix3.DWIBiasCorrect(use_ants=True), name='bias_correct')
+
+    def get_b0_mask_fn(b0_file):
+        import os
+        import nibabel as nib
+        from nipype.utils.filemanip import fname_presuffix
+        from dipy.segment.mask import median_otsu
+
+        mask_file = fname_presuffix(
+            b0_file, suffix='_mask', newpath=os.path.abspath('.')
+        )
+        img = nib.load(b0_file)
+        data, aff = img.get_data(), img.affine
+        _, mask = median_otsu(data, 2, 1)
+        nib.Nifti1Image(mask.astype(float), aff).to_filename(mask_file)
+        return mask_file
+
+    b0mask_node = pe.Node(
+        niu.Function(
+            input_names=['b0_file'], output_names=['mask_file'], function=get_b0_mask_fn),
+        name='getB0Mask')
+
+    dwi_tensor_wf = init_dwi_tensor_wf()
+
+    dwi_wf.connect([
+        (inputnode, gen_idx, [('dwi_file', 'in_file')]),
+        (inputnode, acqp, [('dwi_file', 'in_file')]),
+        (inputnode, dwi_artifacts_wf, [('dwi_file', 'inputnode.dwi_file')]),
+        (dwi_artifacts_wf, avg_b0_0, [('outputnode.out_file', 'in_dwi')]),
+        (dwi_artifacts_wf, dwi_eddy_wf, [('outputnode.out_file', 'inputnode.dwi_file')]),
+        (inputnode, avg_b0_0, [('bval_file', 'in_bval')]),
+        (avg_b0_0, bet_dwi0, [('out_file', 'in_file')]),
+        (inputnode, dwi_eddy_wf, [('bval_file', 'inputnode.bval_file'),
+                                  ('bvec_file', 'inputnode.bvec_file')]),
+        (bet_dwi0, dwi_eddy_wf, [('mask_file', 'inputnode.mask_file')]),
+        (gen_idx, dwi_eddy_wf, [('out_file', 'inputnode.index')]),
+        (dwi_eddy_wf, denoise_eddy, [('outputnode.out_file', 'in_file')]),
+        (dwi_eddy_wf, bias_correct, [('outputnode.out_file', 'in_file'),
+                                     ('outputnode.out_bvec', 'in_bvec')]),
+        (inputnode, bias_correct, [('bval_file', 'in_bval')]),
+        (bias_correct, fslroi, [('out_file', 'in_file')]),
+        (fslroi, b0mask_node, [('roi_file', 'b0_file')]),
+        (dwi_eddy_wf, eddy_quad, [(('outputnode.out_file', get_path), 'base_name'),
+                                  (('outputnode.out_file', get_qc_path), 'output_dir'),
+                                  ('outputnode.out_bvec', 'bvec_file')]),
+        (inputnode, eddy_quad, [('bval_file', 'bval_file')]),
+        (b0mask_node, eddy_quad, [('mask_file', 'mask_file')]),
+        (gen_idx, eddy_quad, [('out_file', 'idx_file')]),
+        (bias_correct, dwi_tensor_wf, [('out_file', 'inputnode.dwi_file')]),
+        (b0mask_node, dwi_tensor_wf, [('mask_file', 'inputnode.mask_file')]),
+        (inputnode, dwi_tensor_wf, [('bval_file', 'inputnode.bval_file')]),
+        (dwi_eddy_wf, dwi_tensor_wf, [('outputnode.out_bvec', 'inputnode.bvec_file')])
+    ])
+
+    if acqp_file:
+        eddy_inputspec = dwi_eddy_wf.get_node('inputnode')
+        eddy_inputspec.inputs.acqp = acqp_file
+        eddy_quad.inputs.param_file = acqp_file
+    else:
+        dwi_wf.connect([
+            (acqp, dwi_eddy_wf, [('out_file', 'inputnode.acqp')]),
+            (acqp, eddy_quad, [('out_file', 'param_file')])
+        ])
+
+    # Decide what ecc will take: topup or fmap
+    # If epi files detected
+    if fmaps:
+        if fmap['suffix'] == 'epi':
+            dwi_wf.connect([
+                (dwi_sdc_wf, dwi_eddy_wf, [('outputnode.out_topup', 'inputnode.topup_fieldcoef'),
+                                           ('outputnode.out_movpar', 'inputnode.topup_movpar')])
+            ])
+    # Otherwise (fieldmaps)
+        if any(s in fmap['suffix'] for s in ['fieldmap', 'phasediff', 'phase1']):
+            dwi_wf.connect([
+                (bet_dwi0, dwi_sdc_wf, [('out_file', 'inputnode.dwi_ref_brain')]),
+                (dwi_sdc_wf, dwi_eddy_wf, [('outputnode.out_fmap', 'inputnode.fieldmap_file')])
+            ])
+
+    dwi_wf.connect([
+        (dwi_eddy_wf, outputnode, [('outputnode.out_file', 'out_dwi')]),
+        (inputnode, outputnode, [('bval_file', 'out_bval')]),
+        (dwi_eddy_wf, outputnode, [('outputnode.out_bvec', 'out_bvec')]),
+        (gen_idx, outputnode, [('out_file', 'index')]),
+        (acqp, outputnode, [('out_file', 'acq_params')]),
+        (b0mask_node, outputnode, [('mask_file', 'out_mask')]),
+        (avg_b0_0, outputnode, [('out_file', 'out_b0_pre')]),
+        (bet_dwi0, outputnode, [('mask_file', 'out_b0_mask_pre')]),
+        (eddy_quad, outputnode, [('qc_json', 'out_eddy_quad_json'),
+                                 ('qc_pdf', 'out_eddy_quad_pdf')]),
+        (dwi_tensor_wf, outputnode, [('outputnode.FA_file', 'out_dtifit_FA'),
+                                     ('outputnode.MD_file', 'out_dtifit_MD'),
+                                     ('outputnode.AD_file', 'out_dtifit_AD'),
+                                     ('outputnode.RD_file', 'out_dtifit_RD'),
+                                     ('outputnode.V1_file', 'out_dtifit_V1'),
+                                     ('outputnode.sse_file', 'out_dtifit_sse')])
+    ])
+
+    return dwi_wf
+
+
+def _get_wf_name(dwi_fname):
+    """
+    Derives the workflow name from the supplied dwi file.
+    """
+
+    from nipype.utils.filemanip import split_filename
+
+    fname = split_filename(dwi_fname)[1]
+    fname_nosub = '_'.join(fname.split('_')[1:])
+    name = 'dwi_preproc_' + fname_nosub.replace('.', '_').replace(' ', '').replace(
+        '-', '_'
+    ).replace('_dwi', '_wf')
+
+    return name
diff --git a/dmriprep/workflows/dwi/eddy.py b/dmriprep/workflows/dwi/eddy.py
new file mode 100644
index 00000000..647edfe7
--- /dev/null
+++ b/dmriprep/workflows/dwi/eddy.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+
+"""
+Head motion, eddy current distortion and EPI distortion correction
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_dwi_eddy_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import fsl, utility as niu
+from numba import cuda
+
+
+def init_dwi_eddy_wf(omp_nthreads, sdc_method=None):
+    """
+    This workflow runs eddy on the input dwi image.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from dmriprep.workflows.dwi import init_dwi_eddy_wf
+        wf = init_dwi_eddy_wf(omp_nthreads=1, sdc_method=None)
+
+    **Parameters**
+
+        omp_nthreads: int
+            Number of threads to run eddy
+        sdc_type: str
+            Synthetic distortion correction method (may include 'fieldmap, 'topup')
+
+    **Inputs**
+
+        dwi_file
+            dwi NIfTI file
+        bvec_file
+            bvec file
+        bval_file
+            bval file
+        mask_file
+            brain mask file
+        fieldmap_file
+            fieldmap file
+        topup_fieldcoef
+            topup file containing field coefficients
+        topup_movpar
+            topup movpar.txt file
+        acqp
+            acquisition parameters file
+        index
+            index file
+
+    **Outputs**
+
+        out_file
+            output eddy-corrected dwi image
+        out_bvec
+            output rotated bvecs after eddy correction
+
+    """
+
+    wf = pe.Workflow(name='dwi_eddy_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=['dwi_file',
+                                                      'bvec_file',
+                                                      'bval_file',
+                                                      'mask_file',
+                                                      'fieldmap_file',
+                                                      'topup_fieldcoef',
+                                                      'topup_movpar',
+                                                      'acqp',
+                                                      'index']),
+                       name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_bvec']),
+                         name='outputnode')
+
+    ecc = pe.Node(
+        fsl.Eddy(num_threads=omp_nthreads, repol=True, cnr_maps=True, residuals=True),
+        name='fsl_eddy')
+
+    try:
+        if cuda.gpus:
+            ecc.inputs.use_cuda = True
+    except:
+        ecc.inputs.use_cuda = False
+
+    wf.connect([
+        (inputnode, ecc, [('dwi_file', 'in_file'),
+                          ('bval_file', 'in_bval'),
+                          ('bvec_file', 'in_bvec'),
+                          ('acqp', 'in_acqp'),
+                          ('index', 'in_index'),
+                          ('mask_file', 'in_mask')]),
+        (ecc, outputnode, [('out_corrected', 'out_file'),
+                           ('out_rotated_bvecs', 'out_bvec')])
+    ])
+
+    if sdc_method == 'fieldmap':
+        wf.connect([
+            (inputnode, ecc, [('fieldmap_file', 'field')])
+        ])
+
+    if sdc_method == 'topup':
+        wf.connect([
+            (inputnode, ecc, [('topup_fieldcoef', 'in_topup_fieldcoef'),
+                              ('topup_movpar', 'in_topup_movpar')])
+        ])
+
+    return wf
diff --git a/dmriprep/workflows/dwi/outputs.py b/dmriprep/workflows/dwi/outputs.py
new file mode 100644
index 00000000..2a755fb5
--- /dev/null
+++ b/dmriprep/workflows/dwi/outputs.py
@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+
+"""
+Output dwi
+^^^^^^^^^^
+
+.. autofunction:: init_dwi_derivatives_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import io as nio, utility as niu
+
+
+def init_dwi_derivatives_wf(subject_id, session_id, output_folder):
+
+    output_wf = pe.Workflow(name="output_wf")
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                "subject_id",
+                "session_id",
+                "dwi",
+                "bval",
+                "bvec",
+                "index",
+                "acq_params",
+                "mask",
+                "b0",
+                "b0_mask",
+                "eddy_quad_json",
+                "eddy_quad_pdf",
+                "dtifit_FA",
+                "dtifit_MD",
+                "dtifit_AD",
+                "dtifit_RD",
+                "dtifit_V1",
+                "dtifit_sse",
+                "noise",
+            ]
+        ),
+        name="inputnode",
+    )
+
+    def build_path(output_folder, subject_id, session_id):
+        import os
+
+        return os.path.join(
+            output_folder, "dmriprep", "sub-" + subject_id, "ses-" + session_id, "dwi"
+        )
+
+    concat = pe.Node(
+        niu.Function(
+            input_names=["output_folder", "subject_id", "session_id"],
+            output_names=["built_folder"],
+            function=build_path,
+        ),
+        name="build_path",
+    )
+
+    datasink = pe.Node(nio.DataSink(), name="datasink")
+
+    output_wf.connect(
+        [
+            (
+                inputnode,
+                concat,
+                [
+                    ("subject_id", "subject_id"),
+                    ("session_id", "session_id"),
+                    ("output_folder", "output_folder"),
+                ],
+            ),
+            (concat, datasink, [("built_folder", "base_directory")]),
+            (
+                inputnode,
+                datasink,
+                [
+                    ("dwi", "@result.@dwi"),
+                    ("bvec", "@result.@bvec"),
+                    ("bval", "@result.@bval"),
+                    ("mask", "@result.@mask"),
+                    ("b0", "@result.@b0"),
+                    ("eddy_quad_json", "@result.@eddy_quad_json"),
+                    ("eddy_quad_pdf", "@result.@eddy_quad_pdf"),
+                    ("dtifit_FA", "@result.@fa"),
+                    ("dtifit_MD", "@result.@md"),
+                    ("dtifit_AD", "@result.@ad"),
+                    ("dtifit_RD", "@result.@rd"),
+                    ("dtifit_V1", "@result.@v1"),
+                    ("dtifit_sse", "@result.@sse"),
+                ],
+            ),
+        ]
+    )
+
+    return output_wf
diff --git a/dmriprep/workflows/dwi/tensor.py b/dmriprep/workflows/dwi/tensor.py
new file mode 100644
index 00000000..92bea601
--- /dev/null
+++ b/dmriprep/workflows/dwi/tensor.py
@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+
+"""
+Tensor fitting
+^^^^^^^^^^^^^^
+
+.. autofunction:: init_dwi_tensor_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import fsl, utility as niu
+
+
+def init_dwi_tensor_wf():
+    """
+    This workflow fits a tensor to the input dwi image.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from dmriprep.workflows.dwi import init_dwi_tensor_wf
+        wf = init_dwi_tensor_wf()
+
+    **Inputs**
+
+        dwi_file
+            dwi NIfTI file
+        bvec_file
+            bvec file
+        bval_file
+            bval file
+        mask_file
+            brain mask file
+
+    **Outputs**
+
+        FA_file
+            FA
+        MD_file
+            MD
+        AD_file
+            AD
+        RD_file
+            RD
+        V1_file
+            V1
+        sse_file
+            diffusion tensor fit residuals
+
+    """
+
+    wf = pe.Workflow(name='dwi_tensor_wf')
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=['dwi_file', 'bvec_file', 'bval_file', 'mask_file']),
+        name='inputnode')
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=['FA_file', 'MD_file', 'AD_file', 'RD_file', 'V1_file', 'sse_file']),
+        name='outputnode')
+
+    dtifit = pe.Node(fsl.DTIFit(save_tensor=True, sse=True), name='dtifit')
+
+    add_l2_l3 = pe.Node(fsl.BinaryMaths(operation='add'), name='add_l2_l3')
+
+    calc_rd = pe.Node(fsl.BinaryMaths(operand_value=2, operation='div'), name='calc_rd')
+
+    wf.connect([
+        (inputnode, dtifit, [('dwi_file', 'dwi'),
+                             ('bvec_file', 'bvecs'),
+                             ('bval_file', 'bvals'),
+                             ('mask_file', 'mask')]),
+        (dtifit, add_l2_l3, [('L2', 'in_file'),
+                             ('L3', 'operand_file')]),
+        (add_l2_l3, calc_rd, [('out_file', 'in_file')]),
+        (calc_rd, outputnode, [('out_file', 'RD_file')]),
+        (dtifit, outputnode, [('FA', 'FA_file'),
+                              ('MD', 'MD_file'),
+                              ('L1', 'AD_file'),
+                              ('V1', 'V1_file'),
+                              ('sse', 'sse_file')])
+    ])
+
+    return wf
diff --git a/dmriprep/workflows/dwi/util.py b/dmriprep/workflows/dwi/util.py
new file mode 100644
index 00000000..0ebf58c7
--- /dev/null
+++ b/dmriprep/workflows/dwi/util.py
@@ -0,0 +1,179 @@
+# -*- coding: utf-8 -*-
+
+"""
+Utility workflows
+^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_dwi_concat_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+
+
+def init_dwi_concat_wf(layout):
+    """
+    This workflow concatenates a list of dwi images as well as their associated
+    bvecs and bvals.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from collections import namedtuple
+        from dmriprep.workflows.dwi import init_dwi_concat_wf
+        BIDSLayout = namedtuple('BIDSLayout', ['root'])
+        wf = init_dwi_concat_wf(layout=BIDSLayout('.'))
+
+    **Parameters**
+
+        layout : BIDSLayout object
+            BIDS dataset layout
+
+    **Inputs**
+
+        ref_file
+            reference dwi NIfTI file
+        dwi_list : :obj:`list`
+            list of dwi NIfTI files
+
+    **Outputs**
+
+        dwi_file : :obj:`str`
+            concatenated dwi NIfTI file
+        bvec_file
+            concatenated bvec file
+        bval_file
+            concatenated bval file
+
+    """
+
+    wf = pe.Workflow(name='dwi_concat_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(
+        fields=['ref_file', 'dwi_list']),
+        name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(
+        fields=['dwi_file', 'bvec_file', 'bval_file']),
+        name='outputnode')
+
+    def gather_bvec_bval(layout, dwi_list):
+        bvec_list = [layout.get_bvec(bvec) for bvec in dwi_list]
+        bval_list = [layout.get_bval(bval) for bval in dwi_list]
+        return bvec_list, bval_list
+
+    gather_bvec_bval = pe.Node(
+        niu.Function(
+            input_names=['layout', 'dwi_list'],
+            output_names=['bvec_list', 'bval_list'],
+            function=gather_bvec_bval
+        ),
+        name='gather_bvec_bval')
+    gather_bvec_bval.inputs.layout = layout
+
+    def concat_dwis(ref_file, dwi_list):
+        import os
+        import numpy as np
+        from nipype.utils.filemanip import fname_presuffix
+        import nibabel as nib
+        from nilearn.image import concat_imgs
+
+        out_file = fname_presuffix(
+            ref_file,
+            newpath=os.path.abspath('.')
+        )
+
+        dwi_data = [nib.load(dwi) for dwi in dwi_list]
+
+        new_nii = concat_imgs(dwi_data)
+
+        hdr = dwi_data[0].header.copy()
+        hdr.set_data_shape(new_nii.shape)
+        hdr.set_xyzt_units('mm')
+        hdr.set_data_dtype(np.float32)
+        nib.Nifti1Image(new_nii.get_data(), dwi_data[0].affine, hdr).to_filename(out_file)
+        return out_file
+
+    concat_dwis = pe.Node(
+        niu.Function(
+            input_names=['ref_file', 'dwi_list'],
+            output_names=['out_file'],
+            function=concat_dwis
+        ),
+        name='concat_dwis')
+
+    def concat_bvecs(ref_file, bvec_list):
+        import os
+        import numpy as np
+        from nipype.utils.filemanip import fname_presuffix
+
+        out_file = fname_presuffix(
+            ref_file,
+            suffix='.bvec',
+            newpath=os.path.abspath('.'),
+            use_ext=False
+        )
+
+        bvec_vals = []
+        for bvec in bvec_list:
+            bvec_vals.append(np.genfromtxt(bvec))
+        np.savetxt(out_file,
+                   np.concatenate((bvec_vals), axis=1),
+                   fmt='%.4f',
+                   delimiter=' ')
+        return out_file
+
+    concat_bvecs = pe.Node(
+        niu.Function(
+            input_names=['ref_file', 'bvec_list'],
+            output_names=['out_file'],
+            function=concat_bvecs
+        ),
+        name='concat_bvecs')
+
+    def concat_bvals(ref_file, bval_list):
+        import os
+        import numpy as np
+        from nipype.utils.filemanip import fname_presuffix
+
+        out_file = fname_presuffix(
+            ref_file,
+            suffix='.bval',
+            newpath=os.path.abspath('.'),
+            use_ext=False
+        )
+
+        bval_vals = []
+        for bval in bval_list:
+            bval_vals.append(np.genfromtxt(bval))
+        np.savetxt(out_file,
+                   np.concatenate((bval_vals), axis=0),
+                   fmt='%i',
+                   delimiter=' ',
+                   newline=' ')
+        return out_file
+
+    concat_bvals = pe.Node(
+        niu.Function(
+            input_names=['ref_file', 'bval_list'],
+            output_names=['out_file'],
+            function=concat_bvals
+        ),
+        name='concat_bvals')
+
+    wf.connect([
+        (inputnode, gather_bvec_bval, [('dwi_list', 'dwi_list')]),
+        (inputnode, concat_dwis, [('ref_file', 'ref_file'),
+                                  ('dwi_list', 'dwi_list')]),
+        (inputnode, concat_bvecs, [('ref_file', 'ref_file')]),
+        (gather_bvec_bval, concat_bvecs, [('bvec_list', 'bvec_list')]),
+        (inputnode, concat_bvals, [('ref_file', 'ref_file')]),
+        (gather_bvec_bval, concat_bvals, [('bval_list', 'bval_list')]),
+        (concat_dwis, outputnode, [('out_file', 'dwi_file')]),
+        (concat_bvecs, outputnode, [('out_file', 'bvec_file')]),
+        (concat_bvals, outputnode, [('out_file', 'bval_file')])
+    ])
+
+    return wf
diff --git a/dmriprep/workflows/fieldmap/__init__.py b/dmriprep/workflows/fieldmap/__init__.py
new file mode 100644
index 00000000..78993c0a
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/__init__.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+.. _sdc_estimation :
+Fieldmap estimation and unwarping workflows
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+.. automodule:: dmriprep.workflows.fieldmap.base
+    :members:
+    :undoc-members:
+    :show-inheritance:
+.. automodule:: dmriprep.workflows.fieldmap.fmap
+    :members:
+    :undoc-members:
+    :show-inheritance:
+.. automodule:: dmriprep.workflows.fieldmap.phdiff
+    :members:
+    :undoc-members:
+    :show-inheritance:
+.. automodule:: dmriprep.workflows.fieldmap.pepolar
+    :members:
+    :undoc-members:
+    :show-inheritance:
+.. automodule:: dmriprep.workflows.fieldmap.syn
+    :members:
+    :undoc-members:
+    :show-inheritance:
+"""
+
+from .base import init_sdc_wf
+from .fmap import init_fmap_wf
+#from .unwarp import init_sdc_unwarp_wf, init_fmap_unwarp_report_wf
+from .phdiff import init_phdiff_wf
+from .pepolar import init_pepolar_wf
+from .syn import init_syn_sdc_wf
+
+__all__ = [
+    'init_sdc_wf',
+    'init_pepolar_wf',
+    'init_fmap_wf',
+    'init_phdiff_wf'
+    #'init_sdc_unwarp_wf',
+    #'init_fmap_unwarp_report_wf',
+    'init_syn_sdc_wf'
+]
diff --git a/dmriprep/workflows/fieldmap/base.py b/dmriprep/workflows/fieldmap/base.py
new file mode 100644
index 00000000..040b2881
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/base.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import fsl, utility as niu
+from nipype import logging
+
+from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+
+LOGGER = logging.getLogger('nipype.workflow')
+FMAP_PRIORITY = {
+    'epi': 0,
+    'fieldmap': 1,
+    'phasediff': 2,
+    'phase': 3,
+    'syn': 4
+}
+
+
+def init_sdc_wf(
+    fmaps,
+    metadata,
+    layout,
+    ignore
+):
+
+    sdc_wf = Workflow(name='sdc_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(
+        fields=['dwi_ref', 'dwi_ref_brain', 'dwi_mask',
+                't1_brain', 'std2anat_xfm', 'template', 'templates']),
+        name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(
+        fields=['dwi_ref', 'dwi_ref_brain', 'dwi_mask',
+                'out_fmap', 'out_topup', 'out_movpar', 'out_enc_file',
+                'out_warp', 'syn_dwi_ref', 'method']),
+        name='outputnode')
+
+    fmaps.sort(key=lambda fmap: FMAP_PRIORITY[fmap['suffix']])
+
+    # No fieldmaps
+    if not fmaps or 'fieldmaps' in ignore:
+        sdc_wf.__postdesc__ = """\
+Susceptibility distortion correction (SDC) has been skipped because the
+dataset does not contain extra field map acquisitions correctly described
+with metadata, and the experimental SDC-SyN method was not explicitly selected.
+"""
+
+    outputnode.inputs.method = 'None'
+    sdc_wf.connect([
+        (inputnode, outputnode, [('dwi_ref', 'dwi_ref'),
+                                 ('dwi_ref_brain', 'dwi_ref_brain'),
+                                 ('dwi_mask', 'dwi_mask')])
+    ])
+    return sdc_wf
+
+    sdc_wf.__postdesc__ = """\
+Creating fieldmap
+"""
+
+    fmap = fmaps[0]
+
+    # topup
+    if fmap['suffix'] == 'epi':
+        from .pepolar import init_pepolar_wf
+        outputnode.inputs.method = 'topup'
+
+        epi_fmaps = [
+            (fmap_['epi'], fmap_['metadata']['PhaseEncodingDirection'])
+            for fmap_ in fmaps
+            if fmap_['suffix'] == 'epi'
+        ]
+
+        pepolar_wf = init_pepolar_wf(metadata, epi_fmaps)
+
+        sdc_wf.connect([
+            (inputnode, pepolar_wf, [('dwi_ref_brain', 'inputnode.b0_stripped')]),
+            (pepolar_wf, outputnode, [('outputnode.out_topup', 'out_topup'),
+                                      ('outputnode.out_movpar', 'out_movpar'),
+                                      ('outputnode.out_enc_file', 'out_enc_file'),
+                                      ('outputnode.out_fmap', 'out_fmap')])
+        ])
+
+    # fieldmap
+    elif fmap['suffix'] == 'fieldmap':
+        from .fmap import init_fmap_wf
+        outputnode.inputs.method = 'fieldmap'
+
+        fmap_wf = init_fmap_wf()
+        fmap_wf.inputs.inputnode.fieldmap = fmap['fieldmap']
+        fmap_wf.inputs.inputnode.magnitude = fmap['magnitude']
+
+        mag_flirt = pe.Node(fsl.FLIRT(dof=6), name='magFlirt')
+
+        fmap_flirt = pe.Node(fsl.FLIRT(apply_xfm=True), name='fmapFlirt')
+
+        sdc_wf.connect([
+            (inputnode, mag_flirt, [('dwi_ref_brain', 'reference')]),
+            (fmap_wf, mag_flirt, [('outputnode.fmap_ref', 'in_file')]),
+            (inputnode, fmap_flirt, [('dwi_ref_brain', 'reference')]),
+            (fmap_wf, fmap_flirt, [('outputnode.fmap', 'in_file')]),
+            (mag_flirt, fmap_flirt, [('out_matrix_file', 'in_matrix_file')]),
+            (fmap_flirt, outputnode, [('out_file', 'out_fmap')])
+        ])
+
+    elif fmap['suffix'] in ('phasediff', 'phase'):
+        from .phasediff import init_phase_wf, init_phdiff_wf
+        from .fmap import init_fmap_wf
+        outputnode.inputs.method = 'fieldmap'
+
+        if fmap['suffix'] == 'phasediff':
+            phase_wf = init_phdiff_wf()
+            phase_wf.inputs.inputnode.phasediff = fmap['phasediff']
+            phase_wf.inputs.inputnode.magnitude1 = [
+                fmap_
+                for key, fmap_ in sorted(fmap.items())
+                if key.startswith('magnitude1')
+            ][0]
+            phase_wf.inputs.inputnode.phases_meta = layout.get_metadata(
+                phase_wf.inputs.inputnode.phasediff)
+
+            fmap_wf = init_fmap_wf()
+
+            sdc_wf.connect([
+                (inputnode, fmap_wf, [('dwi_ref_brain', 'inputnode.b0_stripped')]),
+                (phase_wf, fmap_wf, [('outputnode.out_fmap', 'inputnode.fieldmap')]),
+                (phase_wf, fmap_wf, [('outputnode.out_mag', 'inputnode.magnitude')]),
+                (fmap_wf, outputnode, [('outputnode.out_fmap', 'out_fmap')])
+            ])
+
+        elif fmap['suffix'] == 'phase':
+            phase_wf = init_phase_wf()
+            phase_wf.inputs.inputnode.phasediff = [
+                fmap['phase1'],
+                fmap['phase2']
+            ]
+            phase_wf.inputs.inputnode.magnitude1 = [
+                fmap_
+                for key, fmap_ in sorted(fmap.items())
+                if key.startswith('magnitude1')
+            ][0]
+            phase_wf.inputs.inputnode.phases_meta = [
+                layout.get_metadata(i)
+                for i in phase_wf.inputs.inputnode.phasediff
+            ]
+
+            fmap_wf = init_fmap_wf()
+
+            sdc_wf.connect([
+                (inputnode, fmap_wf, [('dwi_ref_brain', 'inputnode.b0_stripped')]),
+                (phase_wf, fmap_wf, [('outputnode.out_fmap', 'inputnode.fieldmap')]),
+                (phase_wf, fmap_wf, [('outputnode.out_mag', 'inputnode.magnitude')]),
+                (fmap_wf, outputnode, [('outputnode.out_fmap', 'out_fmap')])
+            ])
+
+    else:
+        print('No sdc correction')
+    return sdc_wf
diff --git a/dmriprep/workflows/fieldmap/brainsuite.py b/dmriprep/workflows/fieldmap/brainsuite.py
new file mode 100644
index 00000000..f799eafe
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/brainsuite.py
@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+
+"""
+SDC using BrainSuite
+^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_brainsuite_wf
+
+"""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import brainsuite, utility as niu
+
+
+def init_brainsuite_wf():
+
+    wf = pe.Workflow(name='brainsuite_wf')
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'T1w',
+                'T1w_mask',
+                'dwi_file',
+                'dwi_mask',
+                'bvec_file',
+                'bval_file']),
+        name='inputnode')
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['out_dwi']), name='outputnode'
+    )
+
+    bias_corr = pe.Node(brainsuite.Bfc(minBias=0.5, maxBias=1.5), name='bias_corr')
+
+    # nipype spec doesn't have an output definition yet!
+    bdp = pe.Node(brainsuite.BDP(), name='bdp')
+
+    wf.connect([
+        (inputnode, bias_corr, [('T1w', 'inputMRIFile')]),
+        (bias_corr, bdp, [('outputMRIVolume', 'bfcFile'),
+                          ('T1w_mask', 't1Mask'),
+                          ('dwi_file', 'inputDiffusionData'),
+                          ('dwi_mask', 'dwiMask'),
+                          ('bvec_file', 'BVecBValPair')])
+    ])
diff --git a/dmriprep/workflows/fieldmap/fmap.py b/dmriprep/workflows/fieldmap/fmap.py
new file mode 100644
index 00000000..e0ce8f34
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/fmap.py
@@ -0,0 +1,18 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+.. _sdc_direct_b0 :
+Direct B0 mapping sequences
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+When the fieldmap is directly measured with a prescribed sequence (such as
+:abbr:`SE (spiral echo)`), we only need to calculate the corresponding B-Spline
+coefficients to adapt the fieldmap to the TOPUP tool.
+This procedure is described with more detail `here <https://cni.stanford.edu/\
+wiki/GE_Processing#Fieldmaps>`__.
+This corresponds to the section 8.9.3 --fieldmap image (and one magnitude image)--
+of the BIDS specification.
+"""
+from sdcflows.workflows.fmap import init_fmap_wf
+
+__all__ = ['init_fmap_wf']
diff --git a/dmriprep/workflows/fieldmap/pepolar.py b/dmriprep/workflows/fieldmap/pepolar.py
new file mode 100644
index 00000000..42ae01d8
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/pepolar.py
@@ -0,0 +1,156 @@
+#!/usr/bin/env python
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import fsl, utility as niu
+
+
+def init_pepolar_wf(dwi_meta, epi_fmaps):
+
+    dwi_file_pe = dwi_meta['PhaseEncodingDirection']
+
+    file2dir = dict()
+
+    usable_fieldmaps_matching_pe = []
+    usable_fieldmaps_opposite_pe = []
+
+    for fmap, pe_dir in epi_fmaps:
+        if pe_dir == dwi_file_pe:
+            usable_fieldmaps_matching_pe.append(fmap)
+            file2dir[fmap] = pe_dir
+        elif pe_dir[0] == dwi_file_pe[0]:
+            usable_fieldmaps_opposite_pe.append(fmap)
+            file2dir[fmap] = pe_dir
+
+    if not usable_fieldmaps_opposite_pe:
+        raise Exception(
+            'None of the discovered fieldmaps has '
+            'the right phase encoding direction'
+        )
+
+    wf = pe.Workflow(name='pepolar_wf')
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['b0_stripped']), name='inputnode'
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=['out_topup', 'out_movpar', 'out_fmap', 'out_enc_file']
+        ),
+        name='outputnode',
+    )
+
+    topup_wf = init_topup_wf()
+    topup_wf.inputs.inputnode.altepi_file = usable_fieldmaps_opposite_pe[0]
+    wf.add_nodes([inputnode])
+
+    if not usable_fieldmaps_matching_pe:
+        wf.connect(
+            [(inputnode, topup_wf, [('b0_stripped', 'inputnode.epi_file')])]
+        )
+    else:
+        topup_wf.inputs.inputnode.epi_file = usable_fieldmaps_matching_pe[0]
+
+    epi_list = [
+        topup_wf.inputs.inputnode.epi_file,
+        topup_wf.inputs.inputnode.altepi_file,
+    ]
+    dir_map = {
+        'i': 'x',
+        'i-': 'x-',
+        'j': 'y',
+        'j-': 'y-',
+        'k': 'z',
+        'k-': 'z-',
+    }
+    topup_wf.inputs.inputnode.encoding_directions = [
+        dir_map[file2dir[file]] for file in epi_list
+    ]
+
+    wf.connect(
+        [
+            (
+                topup_wf,
+                outputnode,
+                [
+                    ('outputnode.out_fmap', 'out_fmap'),
+                    ('outputnode.out_movpar', 'out_movpar'),
+                    ('outputnode.out_base', 'out_topup'),
+                    ('outputnode.out_enc_file', 'out_enc_file'),
+                ],
+            )
+        ]
+    )
+
+    return wf
+
+
+def init_topup_wf(output_resolution, acqp_file):
+    from ...interfaces import mrtrix3
+
+    wf = pe.Workflow(name='topup_wf')
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'epi_file',
+                'altepi_file',
+                'encoding_directions',
+                'topup_name',
+                'acqp',
+            ]
+        ),
+        name='inputnode',
+    )
+    inputnode.inputs.topup_name = 'topup_base'
+
+    epi_flirt = pe.Node(fsl.FLIRT(), name='epi_flirt')
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=['out_fmap', 'out_movpar', 'out_base', 'out_enc_file']
+        ),
+        name='outputnode',
+    )
+
+    list_merge = pe.Node(niu.Merge(numinputs=2), name='list_merge')
+
+    merge = pe.Node(fsl.Merge(dimension='t'), name='mergeAPPA')
+
+    # Resize (make optional)
+    resize = pe.Node(mrtrix3.MRResize(), name='epi_resize')
+
+    topup = pe.Node(fsl.TOPUP(), name='topup')
+
+    if acqp_file:
+        wf.connect([(inputnode, topup, [('acqp', 'encoding_file')])])
+    else:
+        topup.inputs.readout_times = [0.05, 0.05]
+        wf.connect([
+            (inputnode, topup, [('encoding_directions', 'encoding_direction')])
+        ])
+
+    wf.connect([
+        (inputnode, list_merge, [('epi_file', 'in1')]),
+        (inputnode, epi_flirt, [('altepi_file', 'in_file'),
+                                ('epi_file', 'reference')]),
+        (epi_flirt, list_merge, [('out_file', 'in2')]),
+        (list_merge, merge, [('out', 'in_files')]),
+        (topup, outputnode, [('out_field', 'out_fmap'),
+                             ('out_movpar', 'out_movpar'),
+                             ('out_enc_file', 'out_enc_file')])
+    ])
+
+    if output_resolution:
+        resize.inputs.voxel_size = output_resolution
+
+        wf.connect([
+            (merge, resize, [('merged_file', 'in_file')]),
+            (resize, topup, [('out_file', 'in_file')])
+        ])
+    else:
+        wf.connect([
+            (merge, topup, [('merged_file', 'in_file')])
+        ])
+
+    return wf
diff --git a/dmriprep/workflows/fieldmap/phdiff.py b/dmriprep/workflows/fieldmap/phdiff.py
new file mode 100644
index 00000000..7bda305c
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/phdiff.py
@@ -0,0 +1,15 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+.. _sdc_phasediff :
+Phase-difference B0 estimation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+The field inhomogeneity inside the scanner (fieldmap) is proportional to the
+phase drift between two subsequent :abbr:`GRE (gradient recall echo)`
+sequence.
+Fieldmap preprocessing workflow for fieldmap data structure
+8.9.1 in BIDS 1.0.0: one phase diff and at least one magnitude image
+"""
+from sdcflows.workflows.phdiff import init_phdiff_wf
+
+__all__ = ['init_phdiff_wf']
diff --git a/dmriprep/workflows/fieldmap/syn.py b/dmriprep/workflows/fieldmap/syn.py
new file mode 100644
index 00000000..606b9faf
--- /dev/null
+++ b/dmriprep/workflows/fieldmap/syn.py
@@ -0,0 +1,21 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+.. _sdc_fieldmapless :
+Fieldmap-less estimation (experimental)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+In the absence of direct measurements of fieldmap data, we provide an (experimental)
+option to estimate the susceptibility distortion based on the ANTs symmetric
+normalization (SyN) technique.
+This feature may be enabled, using the ``--use-syn-sdc`` flag, and will only be
+applied if fieldmaps are unavailable.
+During the evaluation phase, the ``--force-syn`` flag will cause this estimation to
+be performed *in addition to* fieldmap-based estimation, to permit the direct
+comparison of the results of each technique.
+Note that, even if ``--force-syn`` is given, the functional outputs of FMRIPREP will
+be corrected using the fieldmap-based estimates.
+Feedback will be enthusiastically received.
+"""
+from sdcflows.workflows.syn import init_syn_sdc_wf
+
+__all__ = ['init_syn_sdc_wf']
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 00000000..9b90b864
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = python -msphinx
+SPHINXPROJ    = dmriprep
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/docs/changes.rst b/docs/changes.rst
new file mode 100644
index 00000000..1a4f845b
--- /dev/null
+++ b/docs/changes.rst
@@ -0,0 +1,5 @@
+----------
+What's new
+----------
+
+.. include:: ../CHANGES.rst
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100755
index 00000000..463b3576
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# dmriprep documentation build configuration file, created by
+# sphinx-quickstart on Fri Jun  9 13:47:02 2017.
+#
+# This file is execfile()d with the current directory set to its
+# containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+# If extensions (or modules to document with autodoc) are in another
+# directory, add these directories to sys.path here. If the directory is
+# relative to the documentation root, use os.path.abspath to make it
+# absolute, like shown here.
+#
+import os
+import sys
+from datetime import datetime
+
+sys.path.insert(0, os.path.abspath('..'))
+
+# -- General configuration ---------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.doctest',
+    'sphinx.ext.intersphinx',
+    'sphinx.ext.viewcode',
+    'sphinxarg.ext',  # argparse extension
+    'nipype.sphinxext.plot_workflow'
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = 'dmriprep'
+author = 'The dMRIPrep developers'
+copyright = '2018-%s, %s' % (datetime.now().year, author)
+
+
+# The version info for the project you're documenting, acts as replacement
+# for |version| and |release|, also used in various other places throughout
+# the built documents.
+#
+# The short X.Y version.
+version = u'version'
+# The full version, including alpha/beta/rc tags.
+release = u'version'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set 'language' from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This patterns also effect to html_static_path and html_extra_path
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# If true, `todo` and `todoList` produce output, else they produce nothing.
+todo_include_todos = False
+
+
+# -- Options for HTML output -------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'sphinx_rtd_theme'
+
+# Theme options are theme-specific and customize the look and feel of a
+# theme further.  For a list of options available for each theme, see the
+# documentation.
+#
+# html_theme_options = {}
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named 'default.css' will overwrite the builtin 'default.css'.
+html_static_path = ['_static']
+
+
+# -- Options for HTMLHelp output ---------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'dmriprepdoc'
+
+
+# -- Options for LaTeX output ------------------------------------------
+
+latex_elements = {
+    # The paper size ('letterpaper' or 'a4paper').
+    #
+    # 'papersize': 'letterpaper',
+    # The font size ('10pt', '11pt' or '12pt').
+    #
+    # 'pointsize': '10pt',
+    # Additional stuff for the LaTeX preamble.
+    #
+    # 'preamble': '',
+    # Latex figure (float) alignment
+    #
+    # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title, author, documentclass
+# [howto, manual, or own class]).
+latex_documents = [
+    (master_doc, 'dmriprep.tex', 'dMRIPrep Documentation', author, 'manual')
+]
+
+
+# -- Options for manual page output ------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [(master_doc, 'dmriprep', u'dMRIPrep Documentation', [author], 1)]
+
+
+# -- Options for Texinfo output ----------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+#  dir menu entry, description, category)
+texinfo_documents = [
+    (
+        master_doc,
+        'dmriprep',
+        'dMRIPrep Documentation',
+        author,
+        'dmriprep',
+        'One line description of project.',
+        'Miscellaneous',
+    )
+]
diff --git a/docs/dmriprep.interfaces.rst b/docs/dmriprep.interfaces.rst
new file mode 100644
index 00000000..70866b7d
--- /dev/null
+++ b/docs/dmriprep.interfaces.rst
@@ -0,0 +1,30 @@
+dmriprep.interfaces package
+===========================
+
+Submodules
+----------
+
+dmriprep.interfaces.fmap module
+-------------------------------
+
+.. automodule:: dmriprep.interfaces.fmap
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.interfaces.reportlets module
+-------------------------------------
+
+.. automodule:: dmriprep.interfaces.reportlets
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.interfaces
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.rst b/docs/dmriprep.rst
new file mode 100644
index 00000000..459f093d
--- /dev/null
+++ b/docs/dmriprep.rst
@@ -0,0 +1,48 @@
+dmriprep package
+================
+
+Subpackages
+-----------
+
+.. toctree::
+
+   dmriprep.interfaces
+   dmriprep.utils
+   dmriprep.viz
+   dmriprep.workflows
+
+Submodules
+----------
+
+dmriprep.cli module
+-------------------
+
+.. automodule:: dmriprep.cli
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.due module
+-------------------
+
+.. automodule:: dmriprep.due
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.qc module
+------------------
+
+.. automodule:: dmriprep.qc
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.utils.rst b/docs/dmriprep.utils.rst
new file mode 100644
index 00000000..d7abef5a
--- /dev/null
+++ b/docs/dmriprep.utils.rst
@@ -0,0 +1,30 @@
+dmriprep.utils package
+======================
+
+Submodules
+----------
+
+dmriprep.utils.bids module
+--------------------------
+
+.. automodule:: dmriprep.utils.bids
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.utils.hemisphere module
+--------------------------------
+
+.. automodule:: dmriprep.utils.hemisphere
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.utils
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.viz.rst b/docs/dmriprep.viz.rst
new file mode 100644
index 00000000..d972ce7a
--- /dev/null
+++ b/docs/dmriprep.viz.rst
@@ -0,0 +1,22 @@
+dmriprep.viz package
+====================
+
+Submodules
+----------
+
+dmriprep.viz.utils module
+-------------------------
+
+.. automodule:: dmriprep.viz.utils
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.viz
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.workflows.dwi.rst b/docs/dmriprep.workflows.dwi.rst
new file mode 100644
index 00000000..015977ea
--- /dev/null
+++ b/docs/dmriprep.workflows.dwi.rst
@@ -0,0 +1,70 @@
+dmriprep.workflows.dwi package
+==============================
+
+Submodules
+----------
+
+dmriprep.workflows.dwi.artifacts module
+---------------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.artifacts
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.base module
+----------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.base
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.eddy module
+----------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.eddy
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.outputs module
+-------------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.outputs
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.qc module
+--------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.qc
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.tensor module
+------------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.tensor
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.dwi.util module
+----------------------------------
+
+.. automodule:: dmriprep.workflows.dwi.util
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.workflows.dwi
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.workflows.fieldmap.rst b/docs/dmriprep.workflows.fieldmap.rst
new file mode 100644
index 00000000..cd4eda34
--- /dev/null
+++ b/docs/dmriprep.workflows.fieldmap.rst
@@ -0,0 +1,70 @@
+dmriprep.workflows.fieldmap package
+===================================
+
+Submodules
+----------
+
+dmriprep.workflows.fieldmap.ants module
+---------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.ants
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.base module
+---------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.base
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.brainsuite module
+---------------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.brainsuite
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.fmap module
+---------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.fmap
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.pepolar module
+------------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.pepolar
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.phasediff module
+--------------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.phasediff
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.fieldmap.synb0 module
+----------------------------------------
+
+.. automodule:: dmriprep.workflows.fieldmap.synb0
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.workflows.fieldmap
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/dmriprep.workflows.rst b/docs/dmriprep.workflows.rst
new file mode 100644
index 00000000..c2d91d7e
--- /dev/null
+++ b/docs/dmriprep.workflows.rst
@@ -0,0 +1,38 @@
+dmriprep.workflows package
+==========================
+
+Subpackages
+-----------
+
+.. toctree::
+
+   dmriprep.workflows.dwi
+   dmriprep.workflows.fieldmap
+
+Submodules
+----------
+
+dmriprep.workflows.anatomical module
+------------------------------------
+
+.. automodule:: dmriprep.workflows.anatomical
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+dmriprep.workflows.base module
+------------------------------
+
+.. automodule:: dmriprep.workflows.base
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: dmriprep.workflows
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/index.rst b/docs/index.rst
new file mode 100644
index 00000000..d35d1723
--- /dev/null
+++ b/docs/index.rst
@@ -0,0 +1,14 @@
+.. include:: links.rst
+.. include:: ../README.rst
+.. include:: license.rst
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents
+
+   installation
+   changes
+   usage
+   workflows
+   outputs
+   contributing
diff --git a/docs/installation.rst b/docs/installation.rst
new file mode 100644
index 00000000..1dbf0078
--- /dev/null
+++ b/docs/installation.rst
@@ -0,0 +1,89 @@
+.. highlight:: shell
+
+============
+Installation
+============
+
+Docker Container
+================
+
+.. code-block:: console
+
+    $ git clone https://github.com/nipreps/dmriprep
+    $ cd dmriprep
+    $ make docker
+
+If you don't want to log into the docker image:
+
+.. code-block:: console
+
+    $ docker run -ti -v $BIDS_INPUT_DIR:/inputs -v $OUTPUT_DIR:/outputs dmriprep:prod dmriprep /inputs /outputs
+
+If you want to log into the image:
+
+.. code-block:: console
+
+    $ docker run -ti -v $BIDS_INPUT_DIR:/inputs -v $OUTPUT_DIR:/outputs dmriprep:prod
+
+Run this inside the docker image:
+
+.. code-block:: console
+
+    $ dmriprep /inputs /output --participant-label 01
+
+Singularity Container
+=====================
+
+Preparing a Singularity image (Singularity version < 2.5)
+---------------------------------------------------------
+
+Running a Singularity Image
+---------------------------
+
+Manually Prepared Environment (Python 3.5+)
+===========================================
+
+To install dmriprep, run this command in your terminal:
+
+.. code-block:: console
+
+    $ pip install dmriprep
+
+This is the preferred method to install dmriprep, as it will always install the most recent stable release.
+
+If you don't have `pip`_ installed, this `Python installation guide`_ can guide
+you through the process.
+
+.. _pip: https://pip.pypa.io
+.. _Python installation guide: http://docs.python-guide.org/en/latest/starting/installation/
+
+
+From sources
+------------
+
+The sources for dmriprep can be downloaded from the `Github repo`_.
+
+You can either clone the public repository:
+
+.. code-block:: console
+
+    $ git clone git://github.com/nipreps/dmriprep
+
+Or download the `tarball`_:
+
+.. code-block:: console
+
+    $ curl  -OL https://github.com/nipreps/dmriprep/tarball/master
+
+Once you have a copy of the source, you can install it with:
+
+.. code-block:: console
+
+    $ python setup.py install
+
+
+.. _Github repo: https://github.com/nipreps/dmriprep
+.. _tarball: https://github.com/nipreps/dmriprep/tarball/master
+
+External Dependencies
+---------------------
diff --git a/docs/license.rst b/docs/license.rst
new file mode 100644
index 00000000..d48d8a1b
--- /dev/null
+++ b/docs/license.rst
@@ -0,0 +1,11 @@
+License information
+-------------------
+
+We use the Apache license; the full license may be found in the
+`LICENSE <https://github.com/nipreps/dmriprep/blob/master/LICENSE>`_ file
+in the ``dmriprep`` distribution.
+
+All trademarks referenced herein are property of their respective holders.
+
+Copyright (c) 2018-2019, the dMRIPrep developers.
+All rights reserved.
diff --git a/docs/links.rst b/docs/links.rst
new file mode 100644
index 00000000..32dca760
--- /dev/null
+++ b/docs/links.rst
@@ -0,0 +1,6 @@
+.. _Nipype: http://nipype.readthedocs.io/en/latest/
+.. _BIDS: http://bids.neuroimaging.io
+.. _Usage: usage.html
+.. _Installation: installation.html
+.. _workflows: workflows.html
+.. _guidelines: contributing.html
diff --git a/docs/modules.rst b/docs/modules.rst
new file mode 100644
index 00000000..05dd5a95
--- /dev/null
+++ b/docs/modules.rst
@@ -0,0 +1,7 @@
+dmriprep
+========
+
+.. toctree::
+   :maxdepth: 4
+
+   dmriprep
diff --git a/docs/outputs.rst b/docs/outputs.rst
new file mode 100644
index 00000000..1b449523
--- /dev/null
+++ b/docs/outputs.rst
@@ -0,0 +1,46 @@
+Outputs of dMRIPrep
+===================
+
+dMRIPrep generates:
+
+  1. **Visual QA (quality assurance) reports**
+
+  2. **Pre-processed imaging data**
+
+Visual Reports
+--------------
+
+Preprocessed data (dMRIPrep *derivatives*)
+------------------------------------------
+
+Anatomical derivatives are placed in each subject's ``anat`` subfolder:
+
+Diffusion derivatives are stored in the ``dwi`` subfolder.
+
+.. code-block:: console
+
+    dmriprep/
+      sub-001/
+        ses-01/
+          dwi/
+            # preprocessed diffusion weighted images
+            sub-001_ses-01_desc-preproc.dwi.nii.gz
+            sub-001_ses-01_desc-preproc.dwi.bval
+            sub-001_ses-01_desc-preproc.dwi.bvec
+            sub-001_ses-01_desc-preproc.dwi.json
+            # mask
+            sub-001_ses-01_desc-brain_mask.nii.gz
+            # model-derived maps
+            sub-001_ses-01_model-DTI_desc-preproc_FA.nii.gz
+            sub-001_ses-01_model-DTI_desc-preproc_MD.nii.gz
+            sub-001_ses-01_model-DTI_desc-preproc_AD.nii.gz
+            sub-001_ses-01_model-DTI_desc-preproc_RD.nii.gz
+            # directionally-encoded colour (DEC) maps
+            sub-001_ses-01_model-DTI_desc-DEC_FA.nii.gz
+            # grey/white matter masks
+            sub-001_ses-01_desc-aparcaseg_dseg.nii.gz
+            sub-001_ses-01_desc-aseg_dseg.nii.gz
+          anat/
+      dataset_description.json
+      sub-001.html
+      index.html
diff --git a/docs/usage.rst b/docs/usage.rst
new file mode 100644
index 00000000..d4c3e902
--- /dev/null
+++ b/docs/usage.rst
@@ -0,0 +1,60 @@
+.. include:: links.rst
+
+Usage
+-----
+
+Execution and the BIDS format
+=============================
+
+The ``dmriprep`` workflow takes as principal input the path of the dataset
+that is to be processed.
+The input dataset is required to be in valid :abbr:`BIDS (Brain Imaging Data
+Structure)` format, and it must include at least one T1w structural image and
+a dwi series.
+
+.. code-block:: console
+
+    bids
+    └── sub-01
+      ├── anat
+      │   └── sub-01_T1w.json
+      │   └── sub-01_T1w.nii.gz
+      ├── dwi
+      │   ├── sub-01_dwi.bval
+      │   ├── sub-01_dwi.bvec
+      │   ├── sub-01_dwi.json
+      │   └── sub-01_dwi.nii.gz
+      └── fmap
+          ├── sub-01_acq-dwi_dir-AP_epi.json
+          ├── sub-01_acq-dwi_dir-AP_epi.nii.gz
+          ├── sub-01_acq-dwi_dir-PA_epi.json
+          └── sub-01_acq-dwi_dir-PA_epi.nii.gz
+
+We highly recommend that you validate your dataset with the free, online
+`BIDS Validator <http://bids-standard.github.io/bids-validator/>`_.
+
+The exact command to run ``dmriprep`` depends on the Installation_ method.
+The common parts of the command follow the `BIDS-Apps
+<https://github.com/BIDS-Apps>`_ definition.
+
+Example: ::
+
+    dmriprep data/bids out/ participant -w work/
+
+
+Command-Line Arguments
+======================
+
+.. click:: dmriprep.cli:main
+    :prog: dmriprep
+    :show-nested:
+
+Debugging
+=========
+
+Logs and crash files are output into the
+``<output dir>/dmriprep/sub-<participant_label>/log`` directory.
+Information on how to customize and understand these files can be found on the
+`nipype errors and crashes <https://miykael.github.io/nipype_tutorial/notebooks/basic_error_and_crashes.html>`_
+and `nipype debugging <https://miykael.github.io/nipype_tutorial/notebooks/basic_debug.html>`_
+pages.
diff --git a/docs/workflows.rst b/docs/workflows.rst
new file mode 100644
index 00000000..471f5f38
--- /dev/null
+++ b/docs/workflows.rst
@@ -0,0 +1,207 @@
+Processing pipeline details
+===========================
+
+dMRIprep adapts its pipeline depending on what data and metadata are available
+and are used as the input.
+
+A (very) high-level view of the simplest pipeline is presented below.
+
+.. workflow::
+    :graph2use: orig
+    :simple_form: yes
+
+    from collections import namedtuple, OrderedDict
+    from dmriprep.workflows.base import init_single_subject_wf
+    BIDSLayout = namedtuple('BIDSLayout', ['root'])
+    wf = init_single_subject_wf(subject_id='test', session_list=[], name='single_subject_wf', layout=BIDSLayout, output_dir='.', work_dir='.', ignore=['fieldmaps'], b0_thresh=5, output_resolution=(1, 1, 1), bet_dwi=0.3, bet_mag=0.3, omp_nthreads=1, synb0_dir='')
+
+T1w preprocessing
+-----------------
+:mod:`dmriprep.workflows.anatomical.init_anat_preproc_wf`
+
+The T1w preprocessing workflow skull strips the T1w scan for later use in
+tractography or susceptibility distortion correction if using ANTs or
+BrainSuite. It uses ``init_brain_extraction_wf`` from
+`niworkflows <https://github.com/poldracklab/niworkflows/blob/d5904a15201f99ce58a00774d1bfcdd9ea661371/niworkflows/anat/ants.py#L58>`_.
+
+.. workflow::
+    :graph2use: orig
+    :simple_form: yes
+
+    from niworkflows.anat.ants import init_brain_extraction_wf
+    wf = init_brain_extraction_wf()
+
+DWI preprocessing
+-----------------
+:mod:`dmriprep.workflows.dwi.base.init_dwi_preproc_wf`
+
+Preprocessing of DWI files is split into multiple sub-workflows described below.
+
+.. workflow::
+    :graph2use: orig
+    :simple_form: yes
+
+    from collections import namedtuple
+    from dmriprep.workflows.dwi import init_dwi_preproc_wf
+    BIDSLayout = namedtuple('BIDSLayout', ['root'])
+    wf = init_dwi_preproc_wf(layout=BIDSLayout('.'), subject_id='dmripreptest', dwi_file='/madeup/path/sub-01_dwi.nii.gz', metadata={'PhaseEncodingDirection': 'j-', 'TotalReadoutTime': 0.05}, b0_thresh=5, output_resolution=(1, 1, 1), bet_dwi=0.3, bet_mag=0.3, omp_nthreads=1, ignore=['fieldmaps'], synb0_dir='')
+
+Concatenating Scans
+^^^^^^^^^^^^^^^^^^^
+:mod:`dmriprep.workflows.dwi.util.init_dwi_concat_wf`
+
+By default, each scan in the ``dwi`` folder will get preprocessed separately.
+However, there are some cases where multiple scans should be concatenated before
+head motion, eddy current distortion and susceptibility distortion correction
+(eg. single-shell or multi-shell scan acquired in separate runs).
+
+.. code-block:: console
+
+    bids
+    └── sub-01
+      └── dwi
+          ├── sub-01_acq-multishelldir30b1000_dwi.nii.gz
+          ├── sub-01_acq-multishelldir30b3000_dwi.nii.gz
+          ├── sub-01_acq-multishelldir30b4500_dwi.nii.gz
+          └── sub-01_acq-singleshelldir60b1000_dwi.nii.gz
+
+In the above example, the multi-shell scans should be concatenated and the
+single-shell scan should be left alone. This is done by using the
+``--concat_dwis`` flag as shown below:
+
+``--concat_dwis multishelldir30b1000 multishelldir30b3000 multishelldir30b4500``
+
+The unique ``acq`` entities are passed on so this pipeline knows which scans to
+concatenate. The associated bvecs and bvals are also concatenated.
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.dwi import init_dwi_concat_wf
+    wf = init_dwi_concat_wf(ref_file='/madeup/path/sub-01_dwi.nii.gz')
+
+Artifact Removal
+^^^^^^^^^^^^^^^^
+:mod:`dmriprep.workflows.dwi.artifacts.init_dwi_artifacts_wf`
+
+By default, dMRIPrep performs denoising and unringing on each dwi scan. These
+steps can be skipped by passing ``--ignore denoising unringing`` in the command line.
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.dwi import init_dwi_artifacts_wf
+    wf = init_dwi_artifacts_wf(ignore=[], output_resolution=(1, 1, 1))
+
+Motion and Susceptibility Correction
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+:mod:`dmriprep.workflows.dwi.eddy.init_dwi_eddy_wf`
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.dwi import init_dwi_eddy_wf
+    wf = init_dwi_eddy_wf(1, 'topup')
+
+Tensor Estimates
+^^^^^^^^^^^^^^^^
+:mod:`dmriprep.workflows.dwi.tensor.init_dwi_tensor_wf`
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.dwi import init_dwi_tensor_wf
+    wf = init_dwi_tensor_wf()
+
+Susceptibility Distortion Correction (SDC)
+------------------------------------------
+
+Introduction
+^^^^^^^^^^^^
+
+Correction Methods
+^^^^^^^^^^^^^^^^^^
+
+1. topup
+2. fieldmap
+3. phasediff
+4. phase1/phase2
+5. nonlinear registration
+6. synthetic b0
+
+Topup
+~~~~~
+:mod:`dmriprep.workflows.fieldmap.pepolar.init_pepolar_wf`
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.fieldmap import init_pepolar_wf
+    wf = init_pepolar_wf('dmripreptest', 0.3, ['/madeup/path/sub-01_dir-AP_epi.nii.gz', '/madeup/path/sub-01_dir-PA_epi.nii.gz'])
+
+Fieldmap
+~~~~~~~~
+:mod:`dmriprep.workflows.fieldmap.fmap.init_fmap_wf`
+
+All of the fieldmap-related workflows end with ``init_fmap_wf``. Depending on
+the type of fieldmap, other upstream workflows are run. If a phasediff scan is
+given, the phasediff is converted to a proper fieldmap using ``fsl_prepare_fieldmap``.
+If phase1 and phase2 scans are given, they are first converted to a phasediff scan.
+The code for this conversion was derived from `sdcflows <https://github.com/poldracklab/sdcflows/pull/30`_
+and tested on the PNC dataset.
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.fieldmap import init_fmap_wf
+    wf = init_fmap_wf(0.3)
+
+Phasediff
+~~~~~~~~~
+:mod:`dmriprep.workflows.fieldmap.phasediff.init_phdiff_wf`
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.fieldmap import init_phdiff_wf
+    wf = init_phdiff_wf(0.3)
+
+Phase1/Phase2
+~~~~~~~~~~~~~
+:mod:`dmriprep.workflows.fieldmap.phasediff.init_phase_wf`
+
+.. workflow::
+    :graph2use: flat
+    :simple_form: no
+
+    from dmriprep.workflows.fieldmap import init_phase_wf
+    wf = init_phase_wf(0.3)
+
+Synthetic b0
+~~~~~~~~~~~~
+The synb0 method is based off of this `paper <https://www.sciencedirect.com/science/article/abs/pii/S0730725X18306179/>`_. It offers an alternative method of SDC by using deep learning on an anatomical image (T1).
+You can use it in this pipeline by generating the synb0s for the subject(s) and
+passing the bids-like directory containing them to the --synb0_dir parameter.
+To find out how to generate the synb0s, you can visit our
+`forked repo <https://github.com/TIGRLab/Synb0-DISCO>`_.
+Once you have a directory of synb0s (recommended to place as derivatives of
+bids folder, ex. bids/derivatives/synb0/sub-XX), then you are ready to run the
+pipeline using them! Just run dmriprep as you usually would, with ``bids_dir``
+and ``output_dir``, but now add ``--synb0_dir <your_synb0_directory>`` to your command.
+The synb0 acquisition parameters for topup and eddy will be automatically generated in the
+pipeline in the following format:
+
+.. code-block:: console
+
+    0 -1 0 <total_readout_time>
+    0 1 0 0
+
+If you want to overwrite the total_readout_time with one of your own,
+simply add ``--total_readout <new_trt_time>`` to your command.
diff --git a/setup.cfg b/setup.cfg
index f632aae3..27489e72 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -20,15 +20,17 @@ classifiers =
 [options]
 python_requires = >=3.5
 install_requires =
+    dipy ~=1.0.0
     indexed_gzip >=0.8.8
     nibabel >=2.2.1
     nilearn !=0.5.0, !=0.5.1
-    nipype >=1.2.0
+    nipype >=1.2.1
     niworkflows ~= 0.10.1
+    numba ~= 0.44.1
     numpy
     pandas
     psutil >=5.4
-    pybids ~= 0.9.2
+    pybids ~= 0.9.3
     pyyaml
     scikit-image
     smriprep ~= 0.3.0
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 00000000..271e3518
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""Tests for `dmriprep` package."""
+
+import numpy as np
+
+from dmriprep.utils.hemisphere import is_hemispherical
+
+
+def uniform_points_on_sphere(npoints=1, hemisphere=True, rotate=(0, 0, 0)):
+    """Generate random uniform points on a unit (hemi)sphere."""
+    r = 1.0
+    if hemisphere:
+        theta = np.random.uniform(0, np.pi / 2, npoints)
+    else:
+        theta = np.random.uniform(0, np.pi, npoints)
+    phi = np.random.uniform(0, 2 * np.pi, npoints)
+
+    x = r * np.sin(theta) * np.cos(phi)
+    y = r * np.sin(theta) * np.sin(phi)
+    z = r * np.cos(theta)
+
+    vecs = np.stack([x, y, z])
+
+    rot_x = np.array(
+        [
+            [1.0, 0.0, 0.0],
+            [0.0, np.cos(rotate[0]), -np.sin(rotate[0])],
+            [0.0, np.sin(rotate[0]), np.cos(rotate[0])],
+        ]
+    )
+
+    rot_y = np.array(
+        [
+            [np.cos(rotate[1]), 0.0, np.sin(rotate[1])],
+            [0.0, 1.0, 0.0],
+            [-np.sin(rotate[1]), 0.0, np.cos(rotate[1])],
+        ]
+    )
+
+    rot_z = np.array(
+        [
+            [np.cos(rotate[2]), -np.sin(rotate[2]), 0.0],
+            [np.sin(rotate[2]), np.cos(rotate[2]), 0.0],
+            [0.0, 0.0, 1.0],
+        ]
+    )
+
+    vecs = np.dot(rot_z, np.dot(rot_y, np.dot(rot_x, vecs)))
+
+    return vecs.transpose()
+
+
+def test_is_hemispherical():
+    vecs = uniform_points_on_sphere(
+        npoints=100,
+        hemisphere=True,
+        rotate=(
+            np.random.uniform(0, np.pi),
+            np.random.uniform(0, np.pi),
+            np.random.uniform(0, np.pi),
+        ),
+    )
+
+    assert is_hemispherical(vecs)[0]
+    vecs = uniform_points_on_sphere(npoints=100, hemisphere=False)
+    assert not is_hemispherical(vecs)[0]
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 00000000..4347a530
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,19 @@
+[tox]
+envlist = py35, py36, py37, flake8
+
+[travis]
+python =
+    3.7: py37
+    3.6: py36
+    3.5: py35
+
+[testenv:flake8]
+basepython = python
+deps = flake8
+commands = flake8 dmriprep
+
+[testenv]
+setenv =
+    PYTHONPATH = {toxinidir}
+
+commands = python setup.py test
diff --git a/update_changes.sh b/update_changes.sh
new file mode 100644
index 00000000..49c607d5
--- /dev/null
+++ b/update_changes.sh
@@ -0,0 +1,70 @@
+#!/bin/bash
+#
+# Collects the pull-requests since the latest release and
+# arranges them in the CHANGES.rst file.
+#
+# This is a script to be run before releasing a new version.
+# Authored by Oscar Esteban
+#
+# Usage /bin/bash update_changes.sh 1.0.1
+#
+# Copyright (c) 2015-2019, the CRN developers team.
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of dmriprep nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+# Setting      # $ help set
+# set -u         # Treat unset variables as an error when substituting.
+set -x         # Print command traces before executing command.
+
+# Check whether the Upcoming release header is present
+head -1 CHANGES.rst | grep -q Upcoming
+UPCOMING=$?
+if [[ "$UPCOMING" == "0" ]]; then
+    head -n3  CHANGES.rst >> newchanges
+fi
+
+# Elaborate today's release header
+HEADER="$1 ($(date '+%B %d, %Y'))"
+echo $HEADER >> newchanges
+echo $( printf "%${#HEADER}s" | tr " " "=" ) >> newchanges
+echo "" >> newchanges
+
+# Search for PRs since previous release
+git log --grep="Merge pull request" `git describe --tags --abbrev=0`..HEAD --pretty='format:  * %b %s' | sed  's/Merge pull request \#\([^\d]*\)\ from\ .*/(\#\1)/' >> newchanges
+echo "" >> newchanges
+echo "" >> newchanges
+
+# Add back the Upcoming header if it was present
+if [[ "$UPCOMING" == "0" ]]; then
+    tail -n+4 CHANGES.rst >> newchanges
+else
+    cat CHANGES.rst >> newchanges
+fi
+
+# Replace old CHANGES.rst with new file
+mv newchanges CHANGES.rst