--- a
+++ b/docs/source/guide_fitting.rst
@@ -0,0 +1,85 @@
+.. _guide_fitting:
+
+**This guide is still under construction**
+
+Fitting
+-----------
+
+Dosma supports cpu-parallelizable quantitative fitting based on
+`scipy.optimize.curve_fit <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html>`_. 
+
+To perform generic fitting to any array-like object using an arbitrary model function, we can use
+:func:`dosma.curve_fit`. For example, we can use this function to fit an array to a
+monoexponential model ``y = a * exp(b * x)`` using a maximum of 4 workers:
+
+>>> from dosma import curve_fit, monoexponential
+>>> curve_fit(monoexponential, x, y, num_workers=4)
+
+Quantitative fitting is quite common in medical image analysis. For example,
+quantitative MRI (qMRI) has enabled computing voxel-wise relaxation parameter maps
+(e.g. |T2|, |T1rho|, etc.). We can fit a monoexponential model for each voxel across these registered_images,
+where ``tc0`` is the initial guess for parameter :math:`-\frac{1}{b}` in the monoexponential model:
+
+>>> from dosma import MonoExponentialFit
+>>> tc0 = 30.0
+>>> echo_times = np.asarray([10.0, 20.0, 50.0])
+>>> fitter = MonoExponentialFit(tc0=tc0, num_workers=4)
+>>> tc, r2_map = fitter.fit(echo_times, images)
+
+If you don't have a good initial guess for ``tc0`` or expect the initial guess to be dependent on the voxel being fit
+(which is often the case), you can specify that the initial guess should be determined based on results from a
+polynomial fit over the log-linearized form of the monoexponential equation ``log(y) = log(a) - x/tc``:
+
+>>> from dosma import MonoExponentialFit
+>>> tc0 = "polyfit"
+>>> echo_times = np.asarray([10.0, 20.0, 50.0])
+>>> fitter = MonoExponentialFit(tc0=tc0, num_workers=4)
+>>> tc, r2_map = fitter.fit(echo_times, images)
+
+Custom model functions can also be provided and used with :class:`dosma.curve_fit` and :class:`dosma.CurveFitter` (recommended),
+a class wrapper around :class:`dosma.curve_fit` that handles :class:`MedicalVolume` data and supports additional post-processing
+on the fitted parameters. The commands below using :class:`dosma.CurveFitter` and :class:`dosma.curve_fit` are equivalent to the
+``fitter`` above:
+
+>>> from dosma import CurveFitter
+>>> cfitter = CurveFitter(
+... monoexponential, p0=(1.0, -1/tc0), num_workers=4, nan_to_num=0,
+... out_ufuncs=[None, lambda x: -1/x], out_bounds=(0, 100))
+>>> popt, r2_map = cfitter.fit(echo_times, images)
+>>> tc = popt[..., 1]
+
+>>> from dosma import curve_fit
+>>> curve_fit(monoexponential, echo_times, [x.volume for x in images], p0=(1.0, -1/tc0), num_workers=4)
+
+Non-linear curve fitting often requires carefully selected parameter initialization. In cases where
+non-linear curve fitting fails, polynomial fitting may be more effective. Polynomials can be fit to
+the data using :func:`dosma.polyfit` or :class:`dosma.PolyFitter` (recommended),
+which is the polynomial fitting equivalent of ``CurveFitter``. Because polynomial fitting can also be
+done as a single least squares problem, it may also often be faster than standard curve fitting.
+The commands below use ``PolyFitter`` to fit to the log-linearized monoexponential fit
+(i.e. ``log(y) = log(a) + b*x`` to some image data:
+
+>>> from dosma import PolyFitter
+>>> echo_times = np.asarray([10.0, 20.0, 50.0])
+>>> pfitter = PolyFitter(deg=1, nan_to_num=0, out_ufuncs=[None, lambda x: -1/x], out_bounds=(0, 100))
+>>> log_images = [np.log(img) for img in images]
+>>> popt, r2_map = pfitter.fit(echo_times, log_images)
+>>> tc = popt[..., 0]  # note ordering of parameters - see numpy.polyfit for more details.
+
+We can also use the polyfit estimates to initialize the non-linear curve fit. For monoexponential
+fitting, we can do the following:
+
+>>> from dosma import CurveFitter, PolyFitter
+>>> echo_times = np.asarray([10.0, 20.0, 50.0])
+>>> pfitter = PolyFitter(deg=1, r2_threshold=0, num_workers=0)
+>>> log_images = [np.log(img) for img in images]
+>>> popt_pf, _ = pfitter.fit(echo_times, log_images)
+>>> cfitter = CurveFitter(monoexponential, r2_threshold=0.9, nan_to_num=0, out_ufuncs=[None, lambda x: -1/x], out_bounds=(0, 100))
+>>> popt, r2 = cfitter.fit(echo_times, images, p0={"a": popt_pf[..., 1], "b": popt_pf[..., 0]})
+>>> tc = popt[..., 1] 
+
+.. Substitutions
+.. |T2| replace:: T\ :sub:`2`
+.. |T1| replace:: T\ :sub:`1`
+.. |T1rho| replace:: T\ :sub:`1`:math:`{\rho}`
+.. |T2star| replace:: T\ :sub:`2`:sup:`*`