Repository: mehta-lab/recOrder (language: python)

Issue:
[BUG] Rare background correction failure
In using `recOrder` over the last week I had >95% success rate with background corrections. 

Once the scope was running background corrections correctly, it would continue to do so indefinitely. But ~1/10 times I would use an incorrect image as the background correction and give an incorrect reconstruction. I was always able to "fix" this by acquiring and reconstructing again. 

One time I had a dramatic failure:
- I took a background image beside the Kazansky target
- I took a background-corrected acquisition of the Kazansky target (looked as expected)
- I changed samples and took a background image (looked as expected)
- **I moved onto the real sample, took an acquisition, and the result was the Kazansky target**. A repeat acquisition was as expected.

I haven't had time to trace this further, but this is high priority bug. I suspect the recOrder side because the new waveorder doesn't carry any state.
Fix 0.4.0dev tests
This PR against `0.4.0dev` includes the test fixes that should have appeared in #381 (see #383 for CI/CD fix). 


Discussion / hints:

Fix patch:
diff --git a/docs/images/cli_structure.png b/docs/images/cli_structure.png
new file mode 100644
index 00000000..0ddb1dd1
Binary files /dev/null and b/docs/images/cli_structure.png differ
diff --git a/docs/microscope-installation-guide.md b/docs/microscope-installation-guide.md
index d54a00de..b701bd1f 100644
--- a/docs/microscope-installation-guide.md
+++ b/docs/microscope-installation-guide.md
@@ -7,6 +7,22 @@ This guide will walk through a complete recOrder installation consisting of:
 
 Before you start you will need a machine with Windows 10, a Meadowlark DS5020 connected to a liquid crystal device, and a microscope system compatible with `Micromanager`. 
 
+## Install Meadowlark DS5020 and liquid crystals
+
+Start by installing the Meadowlark DS5020 and liquid crystals using the software on the USB stick provided by Meadowlark. You will need to install the USB drivers and CellDrive5000.
+
+**Check your installation versions** by opening CellDrive5000 and double clicking the Meadowlark Optics logo. **We have tested `recOrder == 0.4.0` with "PC software version 1.08" and "Controller firmware version 1.04",** and you will need to upgrade if your software and firmware versions are older. 
+
+To upgrade your "PC software version" use these steps:
+
+- From "Add and remove programs", remove CellDrive5000 and "National Instruments Software".
+- From "Device manager", open the "Meadowlark Optics" group, right click `mlousb`, click "Uninstall device", check "Delete the driver software for this device", and click "Uninstall". Uninstall `Meadowlark Optics D5020 LC Driver` following the same steps.
+- Using the USB stick provided by Meadowlark, reinstall the USB drivers and CellDrive5000. 
+- Confirm that "PC software version" == 1.08
+- **Upgrading users:** you will need to reinstall the Meadowlark device to your micromanager configuration file, because the device driver's name has changed to `MeadowlarkLC`. 
+
+To upgrade your DS5020's firmware, use Meadowlark's "Firmware Updater".
+
 ## Install recOrder software
 
 (Optional but recommended) install [anaconda](https://www.anaconda.com/products/distribution) and create a virtual environment  
@@ -27,11 +43,11 @@ should launch napari (may take 15 seconds on a fresh installation) with the recO
  
 ## Install and configure `Micromanager`
 
-Install `Micromanager 2.0` nightly build `20220920` (https://micro-manager.org/Micro-Manager_Nightly_Builds). 
+Download and install [`Micromanager 2.0` nightly build `20230426` (~150 MB link).](https://download.micro-manager.org/nightly/2.0/Windows/MMSetup_64bit_2.0.1_20230426.exe)
 
-**Note:** We have tested recOrder with `20220920`, but most features will work with newer builds. We recommend testing a minimal installation with `20220920` before testing with a different nightly build or additional device drivers. 
+**Note:** We have tested recOrder with `20230426`, but most features will work with newer builds. We recommend testing a minimal installation with `20230426` before testing with a different nightly build or additional device drivers. 
 
-Before launching `Micromanager`, download the Meadowlark device adapters and calibration files from the [release page](https://github.com/mehta-lab/recOrder/releases/) and place these three unzipped files into your `Micromanager` folder (likely `C:\Program Files\Micro-Manager` or similar). 
+Before launching `Micromanager`, download the USB driver dll from the [release page](https://github.com/mehta-lab/recOrder/releases/) and place this pair of unzipped files into your `Micromanager` folder (likely `C:\Program Files\Micro-Manager` or similar). 
 
 Launch `Micromanager`, open `Devices > Hardware Configuration Wizard...`, and add the `MeadowlarkLcOpenSource` device to your configuration. Confirm your installation by opening `Devices > Device Property Browser...` and confirming that `MeadowlarkLCOpenSource` properties appear. 
 
diff --git a/docs/reconstruction-guide.md b/docs/reconstruction-guide.md
index bbf7990e..cbfd5924 100644
--- a/docs/reconstruction-guide.md
+++ b/docs/reconstruction-guide.md
@@ -1,13 +1,25 @@
 # Automating reconstructions
 
-`recOrder` is undergoing changes to the way that it handles automated reconstructions.
+`recOrder` uses a configuration-file-based command-line interface (CLI) to perform all reconstructions.
 
-`recOrder==0.2.0` had an offline mode with `config.yml` files that were used to configure reconstructions. These config files could be generated and saved via the GUI, and the reconstructions could be run via GUI or CLI. 
+## How can I use `recOrder`'s CLI to perform reconstructions?
+`recOrder`'s CLI is summarized in the following figure:
+<img src="./images/cli_structure.png" align="center">
 
-Although the offline mode had many valuable and convenient features, we found that it had diverged from the online mode and it was difficult to recreate results between online and offline modes. These design limitations led us to the following plan for our upcoming releases:
+The main command `reconstruct` command is composed of two subcommands: `compute-tf` and `apply-inv-tf`. 
 
-`recOrder==0.3.0` (release candidate in the coming weeks) will remove the offline mode and use a set of scripts instead. We recommend modifying the reconstruction scripts in `recOrder/examples/` to recreate and automate reconstructions. 
+A reconstruction can be performed with a single `reconstruct` call. For example:
+```
+recorder reconstruct data.zarr/0/0/0 -c config.yml -o reconstruction.zarr
+```
+Equivalently, a reconstruction can be performed with a `compute-tf` call followed by an `apply-inv-tf` call. For example:
+```
+recorder compute-tf data.zarr/0/0/0 -c config.yml -o tf.zarr
+recorder apply-inv-tf data.zarr/0/0/0 tf.zarr -c config.yml -o reconstruction.zarr
+```
+Computing the transfer function is typically the most expensive part of the reconstruction, so saving a transfer function then applying it to many datasets can save time. 
 
-Although these scripts are not as user friendly as a GUI+CLI solution, we are preparing for a much cleaner solution in `1.0.0`, and we appreciate your patience as we go through this change. We ask users who are having any difficulty with the scripts to [open an issue](https://github.com/mehta-lab/recOrder/issues/new/choose) or [send us an email](mailto:shalin.mehta@czbiohub.org,talon.chandler@czbiohub.org).
+## What types of reconstructions are supported?
+See `/recOrder/examples/` for a list of example configuration files. 
 
-`recOrder==1.0.0` will use a single mode to acquire and reconstruct the data. We are currently planning a refactor that will enable an "acquire once, quickly iterate your reconstruction" workflow, an "acquire now, reconstruct later" workflow, and a "live acquisition, live reconstruction" workflow, among others.  
+TODO: Expand this documentation...need docs for each reconstruction type and parameter.
diff --git a/examples/3D-bf-to-2D-phase.py b/examples/3D-bf-to-2D-phase.py
deleted file mode 100644
index c46b6552..00000000
--- a/examples/3D-bf-to-2D-phase.py
+++ /dev/null
@@ -1,64 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_phase2D,
-)
-from recOrder.compute.phantoms import bf_3D_from_phantom
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-# Option 1: use random data and run this script as is.
-data = bf_3D_from_phantom()  # (Z, Y, X)
-
-# Option 2: load from file
-# reader = read_micromanager("/path/to/ome-tiffs/or/zarr/store/")
-# position, time, channel = 0, 0, 0
-# data = reader.get_array(position)[time, channel, ...]  # read 3D volume
-
-## Set up a reconstructor.
-Z, Y, X = data.shape
-reconstructor_args = {
-    "image_dim": (Y, X),
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "n_slices": Z,  # number of slices in z-stack
-    "z_step_um": 2,  # z-step size in um
-    "in_focus_slice": None,  # integer, set to None to use the central slice
-    "wavelength_nm": 532,
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "NA_illu": 0.2,  # numerical aperture of condenser
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "mode": "2D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-reconstructor = initialize_reconstructor(
-    pipeline="PhaseFromBF", **reconstructor_args
-)
-
-phase2D = reconstruct_phase2D(
-    data, reconstructor, method="Tikhonov", reg_p=1e0
-)
-print(f"Shape of 2D phase data: {np.shape(phase2D)}")
-
-## Save to zarr
-with open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=["Phase"],
-) as dataset:
-    # Write to position "0", with length-one time, channel, and z dimensions
-    dataset["0"] = phase2D[np.newaxis, np.newaxis, np.newaxis, ...]
-
-# These lines open the reconstructed images
-# Alternatively, drag and drop the zarr store into napari and use the recOrder-napari reader.
-v = napari.Viewer()
-v.add_image(data)
-v.add_image(phase2D)
-napari.run()
diff --git a/examples/3D-bf-to-3D-phase.py b/examples/3D-bf-to-3D-phase.py
deleted file mode 100644
index 8ab8b167..00000000
--- a/examples/3D-bf-to-3D-phase.py
+++ /dev/null
@@ -1,64 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_phase3D,
-)
-from recOrder.compute.phantoms import bf_3D_from_phantom
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-# Option 1: use data simulated from a phantom and run this script as is.
-data = bf_3D_from_phantom()  # (Z, Y, X)
-
-# Option 2: load from file
-# reader = read_micromanager('/path/to/ome-tiffs/or/zarr/store/')
-# position, time, channel = 0, 0, 0
-# data = reader.get_array(position)[time, channel, ...]  # read 3D volume
-
-## Set up a reconstructor.
-Z, Y, X = data.shape
-reconstructor_args = {
-    "image_dim": (Y, X),
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "n_slices": Z,  # number of slices in z-stack
-    "z_step_um": 2,  # z-step size in um
-    "wavelength_nm": 532,
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "NA_illu": 0.2,  # numerical aperture of condenser
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-    "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-reconstructor = initialize_reconstructor(
-    pipeline="PhaseFromBF", **reconstructor_args
-)
-
-phase3D = reconstruct_phase3D(
-    data, reconstructor, method="Tikhonov", reg_re=1e-2
-)
-print(f"Shape of 3D phase data: {np.shape(phase3D)}")
-
-## Save to zarr
-with open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=["Phase"],
-) as dataset:
-    # Write to position "0", with length-one time and channel dimensions
-    dataset["0"] = phase3D[np.newaxis, np.newaxis, ...]
-
-# These lines open the reconstructed images
-# Alternatively, drag and drop the zarr store into napari and use the recOrder-napari reader.
-v = napari.Viewer()
-v.add_image(data)
-v.add_image(phase3D)
-napari.run()
diff --git a/examples/3D-fluor-to-3D-density.py b/examples/3D-fluor-to-3D-density.py
deleted file mode 100644
index b4a6665b..00000000
--- a/examples/3D-fluor-to-3D-density.py
+++ /dev/null
@@ -1,61 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_density_from_fluorescence,
-)
-from recOrder.compute.phantoms import fluorescence_from_phantom
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-# Option 1: use data simulated from a phantom and run this script as is.
-data = fluorescence_from_phantom()  # (Z, Y, X)
-
-# Option 2: load from file
-# reader = read_micromanager('/path/to/ome-tiffs/or/zarr/store/')
-# position, time, channel = 0, 0, 0
-# data = reader.get_array(position)[time, channel, ...]  # read 3D volume
-
-## Set up a reconstructor.
-Z, Y, X = data.shape
-reconstructor_args = {
-    "image_dim": (Y, X),
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "n_slices": Z,  # number of slices in z-stack
-    "z_step_um": 2,  # z-step size in um
-    "wavelength_nm": 532,
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-    "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-reconstructor = initialize_reconstructor(
-    pipeline="fluorescence", **reconstructor_args
-)
-
-density = reconstruct_density_from_fluorescence(data, reconstructor, reg=1e-2)
-print(f"Shape of 3D density data: {np.shape(density)}")
-
-## Save to zarr
-with open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=["Density"],
-) as dataset:
-    # Write to position "0", with length-one time and channel dimensions
-    dataset["0"] = density[np.newaxis, np.newaxis, ...]
-
-# These lines open the reconstructed images
-# Alternatively, drag and drop the zarr store into napari and use the recOrder-napari reader.
-v = napari.Viewer()
-v.add_image(data)
-v.add_image(density)
-napari.run()
diff --git a/examples/3D-pol-to-birefringence.py b/examples/3D-pol-to-birefringence.py
deleted file mode 100644
index dfff646d..00000000
--- a/examples/3D-pol-to-birefringence.py
+++ /dev/null
@@ -1,101 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.io.utils import load_bg
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_qlipp_stokes,
-    reconstruct_qlipp_birefringence,
-    reconstruct_phase3D,
-)
-from recOrder.compute.phantoms import pol_3D_from_phantom
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-# Option 1: use random data and run this script as is.
-data, bg_data = pol_3D_from_phantom()  # (C, Z, Y, X) and (C, Y, X)
-C, Z, Y, X = data.shape
-
-# Option 2: load from file
-# reader = read_micromanager('/path/to/ome-tiffs/or/zarr/store/')
-# position, time = 0, 0
-# data = reader.get_array(position)[time, ...]
-# C, Z, Y, X = data.shape
-# bg_data = load_bg("/path/to/recorder/BG", height=Y, width=X)
-
-## Set up a reconstructor.
-reconstructor_args = {
-    "image_dim": (Y, X),
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "n_slices": Z,  # number of slices in z-stack
-    "z_step_um": 2,  # z-step size in um
-    "wavelength_nm": 532,
-    "swing": 0.1,
-    "calibration_scheme": "5-State",  # "4-State" or "5-State"
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "NA_illu": 0.2,  # numerical aperture of condenser
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-    "bg_correction": "local_fit",  # BG correction method: "None", "local_fit", "global"
-    "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-reconstructor = initialize_reconstructor(
-    pipeline="QLIPP", **reconstructor_args
-)
-
-# Reconstruct background Stokes
-bg_stokes = reconstruct_qlipp_stokes(bg_data, reconstructor)
-print(f"Shape of BG Stokes: {np.shape(bg_stokes)}")
-
-# Reconstruct data Stokes w/ background correction
-stokes = reconstruct_qlipp_stokes(data, reconstructor, bg_stokes)
-print(f"Shape of background corrected data Stokes: {np.shape(stokes)}")
-
-# Reconstruct Birefringence from Stokes
-# Shape of the output birefringence will be (C, Z, Y, X) where
-# Channel order = Retardance [nm], Orientation [rad], Brightfield (S0), Degree of Polarization
-birefringence = reconstruct_qlipp_birefringence(stokes, reconstructor)
-birefringence[0] = (
-    birefringence[0] / (2 * np.pi) * reconstructor_args["wavelength_nm"]
-)
-print(f"Shape of birefringence data: {np.shape(birefringence)}")
-
-# Reconstruct Phase3D from S0
-S0 = birefringence[2]
-
-phase3D = reconstruct_phase3D(
-    S0, reconstructor, method="Tikhonov", reg_re=1e-2
-)
-print(f"Shape of 3D phase data: {np.shape(phase3D)}")
-
-## Save to zarr
-with open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=[
-        "Retardance",
-        "Orientation",
-        "BF - computed",
-        "DoP",
-        "Phase",
-    ],
-) as dataset:
-    img = dataset.create_zeros("0", (1, 5, Z, Y, X), dtype=birefringence.dtype)
-    img[0, 0:4] = birefringence
-    img[0, 4] = phase3D
-
-# These lines open the reconstructed images
-# Alternatively, drag and drop the zarr store into napari and use the recOrder-napari reader.
-v = napari.Viewer()
-v.add_image(data)
-v.add_image(phase3D)
-v.add_image(birefringence, contrast_limits=(0, 25))
-v.dims.current_step = (0, 5, 256, 256)
-napari.run()
diff --git a/examples/README.md b/examples/README.md
index 732ef878..adde2a8e 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -1,61 +1,74 @@
-# `recOrder` example scripts
+# `recOrder` CLI examples 
 
-These scripts demonstrate four types of reconstruction:
+`recOrder` uses a configuration-file-based command-line inferface (CLI) to 
+calculate transfer functions and apply these transfer functions to datasets. 
 
-1. 3D brightfield data to 3D phase in `3D-bf-to-3D-phase.py`. 
-2. 3D brightfield data to 2D phase in `3D-bf-to-2D-phase.py`.
-3. 3D + polarization data to retardance, orientation, and 3D phase in `3D-pol-to-birefringence.py`.
-4. 3D fluorescence data to 3D fluorophore density (fluorescence deconvolution) in `3D-fluor-to-3D-density.py`.    
+This page demonstrates `recOrder`'s CLI. 
 
 ## Getting started
 
 ### 1. Check your installation
-First, [install `recOrder`](../docs/software-installation-guide.md) and run any of these scripts with 
+First, [install `recOrder`](../docs/software-installation-guide.md) and run
 ```bash
-python <script-name.py>
+recOrder
+```
+in a shell. If `recOrder` is installed correctly, you will see a usage string and 
+```
+recOrder: Computational Toolkit for Label-Free Imaging
 ```
-Running these scripts without modification will test your installation by running a reconstruction on a simulated phantom. A successful run will open a `napari` window with simulated data and a reconstruction. 
-
-### 2. Reconstruct test data
-Next, [download the test data from zenodo (47 MB)](https://zenodo.org/record/6983916/files/recOrder_test_data.zip?download=1), and modify the script to load the test data. For example, in `3D-bf-to-2D-phase.py`, replace `/path/to/ome-tiffs/or/zarr/store/` with `/path/to/recOrder_test_data/2022_08_04_recOrder_pytest_20x_04NA_BF/2T_3P_16Z_128Y_256X_Kazansky_BF_1`.
-
-### 3. Load and reconstruct your data
 
-Start by loading your data (our readers currently support `.tiff`, `ome.tiff`, `NDTiff`, and `.zarr`, but any `numpy` array will work). Fill in your imaging parameters (or connect your metadata to these parameters), and prototype with this modified script. 
+### 2. Download and convert a test dataset
+Next, [download the test data from zenodo (47 MB)](https://zenodo.org/record/6983916/files/recOrder_test_data.zip?download=1), and convert a dataset to the latest version of `.zarr` with 
+```
+cd /path/to/
+iohub convert -i /path/to/recOrder_test_data/2022_08_04_recOrder_pytest_20x_04NA/2T_3P_16Z_128Y_256X_Kazansky_1/
+-o ./dataset.zarr
+```
 
-We recommend prototyping a reconstruction with a single position and time point so that you can perform initial iterations quickly. 
+You can view the test dataset with
+```
+napari ./dataset.zarr --plugin recOrder-napari
+```
 
-### 4. Sweep your reconstruction parameters
+### 3. Run a reconstruction
+Run an example reconstruction with
+```
+recOrder reconstruct ./dataset.zarr/0/0/0 -c /path/to/recOrder/examples/settings/birefringence-and-phase.yml -o ./reconstruction.zarr
+```
+then view the reconstruction with 
+```
+napari ./reconstruction.zarr --plugin recOrder-napari
+```
 
-You may need to test several parameters to find a value that yields the best results for your application. For example, choosing a regularization parameter is commonly semi-empirical: we recommend choosing a regularization parameter that gives results that aren't too noisy or too smooth. 
+Try modifying the configuration file to see how the regularization parameter changes the results. 
 
-The `sweep-regularization.py` script demonstrates a `3D-bf-to-2D-phase` reconstruction with multiple regularization parameters. We recommend running and understanding this script before modifying your reconstruction script to sweep over regularization or other uncertain parameters to help you settle on a set of reconstruction parameters. 
+### 4. Parallelize over positions or time points
 
-### 5. Reconstruct a multi-modal dataset in a single script
+TODO: @Ed I think we'll use the same strategy as ```mantis deskew```
 
-You may need to perform several reconstructions on a multi-modal dataset. For example, you would like to perform a `3D-fluor-to-3D-density` reconstruction on the fluorescence channels and a `3D-bf-2D-phase` reconstruction on the brightfield channel. 
+## FAQ
+1. **Q: Which configuration file should I use?**
 
-The `multi-modal-recon.py` script demonstrates this type of reconstruction. We recommend running and understanding this script before modifying single-modality reconstructions to run a multi-modal reconstruction. 
+    If you are acquiring:
 
-### 6. Parallelize over positions or time points
+    **3D data with calibrated liquid-crystal polarizers via `recOrder`** use `birefringence.yml`.
 
-Once you've settled on a script that performs a reconstruction, the script can be applied to multiple datasets with a python `for` loop (slowest), `multiprocessing` (faster, see `parallel-reconstruct.py` for an example), or batch processing with an HPC scheduler e.g. `slurm` (fastest). 
+    **3D fluorescence data** use `fluorescence.yml`.
 
-## FAQ
-1. **Q: Which script should I use?**
+    **3D brightfield data** use `phase.yml`.
 
-    If you are acquiring:
+    **Multi-modal data**, start by reconstructing the individual modaliities, each with a single config file and CLI call. Then combine the reconstructions by ***TODO: @Ziwen do can you help me append to the zarrs to help me fix this? ***
 
-    **3D data with calibrated liquid-crystal polarizers via `recOrder`** use `3D-pol-to-birefringence.py`.
+2. **Q: Should I use `reconstruction_dimension` = 2 or 3?
 
-    **3D fluorescence data** use `3D-fluor-to-3D-density.py`.
+    If your downstream processing requires 3D information or if you're unsure, then you should use `reconstruction_dimension = 3`. If your sample is very thin compared to the depth of field of the microscope, if you're in a noise-limited regime, or if your downstream processing requires 2D information, then you should use `reconstruction_dimension = 2`. Empirically, we have found that 2D reconstructions reduce the noise in our reconstructions because it uses 3D information to make a single  estimate for each pixel. 
 
-    **3D brightfield data** use `3D-bf-to-3D-phase.py` or `3D-bf-to-2D-phase.py`, and decide if you need a 3D or 2D phase reconstruction. 
+3. **Q: What regularization parameter should I use?**
 
-    If your downstream processing requires 3D information or if you're unsure, then you should use `3D-bf-to-3D-phase.py`. If your sample is very thin compared to the depth of field of the microscope, if you're in a noise-limited regime, or if your downstream processing requires 2D phase information, then you should use `3D-bf-to-2D-phase.py`. Empirically, we have found that `3D-bf-to-2D-phase.py` reduces the noise in our reconstructions because it uses 3D information to make a single phase estimate for each pixel. 
+    We recommend starting with the defaults then testing over a few orders of magnitude and choosing a result that isn't too noisy or too smooth.
 
-    **Multi-modal data**, start by reconstructing the individual modaliities, then combine the reconstructions using `multi-modal-recon.py` as a guide.  
+### Developers note
 
-2. **Q: What regularization parameter should I use?**
+These configuration files are automatically generated when the tests run. See `/tests/cli_tests/test_settings.py` - `test_generate_example_settings`. 
 
-    We recommend starting with the defaults then testing over a few orders of magnitude and choosing a result that isn't too noisy or too smooth.
+To keep these settings up to date, run `pytest` locally when `cli/settings.py` changes. 
diff --git a/examples/birefringence-and-phase.yml b/examples/birefringence-and-phase.yml
new file mode 100644
index 00000000..fcf30b0c
--- /dev/null
+++ b/examples/birefringence-and-phase.yml
@@ -0,0 +1,30 @@
+input_channel_names:
+- State0
+- State1
+- State2
+- State3
+reconstruction_dimension: 3
+birefringence:
+  transfer_function:
+    swing: 0.1
+  apply_inverse:
+    wavelength_illumination: 0.532
+    background_path: ''
+    remove_estimated_background: false
+    orientation_flip: false
+    orientation_rotate: false
+phase:
+  transfer_function:
+    wavelength_illumination: 0.532
+    yx_pixel_size: 0.325
+    z_pixel_size: 2.0
+    z_padding: 0
+    index_of_refraction_media: 1.3
+    numerical_aperture_detection: 1.2
+    numerical_aperture_illumination: 0.5
+    axial_flip: false
+  apply_inverse:
+    reconstruction_algorithm: Tikhonov
+    regularization_strength: 0.001
+    TV_rho_strength: 0.001
+    TV_iterations: 1
diff --git a/examples/birefringence.yml b/examples/birefringence.yml
new file mode 100644
index 00000000..45ffd268
--- /dev/null
+++ b/examples/birefringence.yml
@@ -0,0 +1,15 @@
+input_channel_names:
+- State0
+- State1
+- State2
+- State3
+reconstruction_dimension: 3
+birefringence:
+  transfer_function:
+    swing: 0.1
+  apply_inverse:
+    wavelength_illumination: 0.532
+    background_path: ''
+    remove_estimated_background: false
+    orientation_flip: false
+    orientation_rotate: false
diff --git a/examples/fluorescence.yml b/examples/fluorescence.yml
new file mode 100644
index 00000000..02d47666
--- /dev/null
+++ b/examples/fluorescence.yml
@@ -0,0 +1,16 @@
+input_channel_names:
+- GFP
+reconstruction_dimension: 3
+fluorescence:
+  transfer_function:
+    yx_pixel_size: 0.325
+    z_pixel_size: 2.0
+    z_padding: 0
+    index_of_refraction_media: 1.3
+    numerical_aperture_detection: 1.2
+    wavelength_emission: 0.507
+  apply_inverse:
+    reconstruction_algorithm: Tikhonov
+    regularization_strength: 0.001
+    TV_rho_strength: 0.001
+    TV_iterations: 1
diff --git a/examples/multi-modal-recon.py b/examples/multi-modal-recon.py
deleted file mode 100644
index c9aadff6..00000000
--- a/examples/multi-modal-recon.py
+++ /dev/null
@@ -1,97 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_phase3D,
-    reconstruct_density_from_fluorescence,
-)
-from recOrder.compute.phantoms import (
-    bf_3D_from_phantom,
-    fluorescence_from_phantom,
-)
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-## Option 1: use data simulated from a phantom and run this script as is.
-bf_data = bf_3D_from_phantom()  # (Z, Y, X)
-fluor_data = fluorescence_from_phantom()  # (Z, Y, X)
-
-## Option 2: load from file. Uncomment all lines with a single #.
-# reader = read_micromanager('/path/to/ome-tiffs/or/zarr/store/')
-# position, time = 0, 0
-
-## Set these by hand or use `reader.channel_names.index('BF')` to read from metadata
-# bf_index = 0
-# fluor_index = 1
-
-## Read
-# bf_data = reader.get_array(position)[time, bf_index, ...]
-# fluor_data = reader.get_array(position)[time, fluor_index, ...]
-
-## Set up the reconstructors.
-
-# Set arguments that are shared by both modalities
-reconstructor_args = {
-    "image_dim": bf_data.shape[1:],  # (Y, X)
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "z_step_um": 2,  # z-step size in um
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "NA_illu": 0.2,  # numerical aperture of condenser
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-    "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-
-# Initialize reconstructors with parameters that are not shared
-bf_reconstructor = initialize_reconstructor(
-    pipeline="PhaseFromBF",
-    wavelength_nm=532,
-    n_slices=bf_data.shape[0],
-    **reconstructor_args
-)
-
-fluor_reconstructor = initialize_reconstructor(
-    pipeline="fluorescence",
-    wavelength_nm=450,
-    n_slices=fluor_data.shape[0],
-    **reconstructor_args
-)
-
-# Reconstruct
-phase = reconstruct_phase3D(
-    bf_data, bf_reconstructor, method="Tikhonov", reg_re=1e-2
-)
-
-density = reconstruct_density_from_fluorescence(
-    fluor_data, fluor_reconstructor, reg=1e-2
-)
-
-## Save to zarr
-with open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=["Phase"],
-) as dataset:
-    # Write to position "0", with length-one time dimension
-    dataset["0"] = phase[np.newaxis, np.newaxis, ...]
-    dataset.append_channel("Density")
-    dataset["0"][0, -1] = density
-
-# These lines open the reconstructed images
-# Alternatively, drag and drop the zarr store into napari and use the recOrder-napari reader.
-v = napari.Viewer()
-v.add_image(bf_data)
-v.add_image(fluor_data)
-v.add_image(phase)
-v.add_image(density)
-v.dims.current_step = (15, 256, 256)
-
-napari.run()
diff --git a/examples/parallel-reconstruct.py b/examples/parallel-reconstruct.py
deleted file mode 100644
index 8cfc4ec3..00000000
--- a/examples/parallel-reconstruct.py
+++ /dev/null
@@ -1,184 +0,0 @@
-import shutil
-import os
-import numpy as np
-import multiprocessing as mp
-from wget import download
-from datetime import datetime
-from iohub import open_ome_zarr
-from iohub.convert import TIFFConverter
-from recOrder.io.utils import load_bg
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_qlipp_stokes,
-    reconstruct_qlipp_birefringence,
-    reconstruct_phase3D,
-)
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-# This example demonstrates a 3D-pol-to-birfringence reconstruction applied
-# to each position and time point of a micromanager dataset. The example
-# uses multiprocessing to apply the reconstructions in parallel.
-
-# This example will download a ~50 MB micro-manager test dataset to your
-# working directory, convert it to a zarr store, then reconstruct physical
-# paramaters from each position and time point.
-
-# Choose the number of processes to parallel over.
-# You should have at least this many cores available with enough RAM
-# available to each process to complete the reconstruction.
-N_processes = 6
-
-
-def precomputation():
-
-    # download test data and convert to zarr
-    # use this section for a test run, then delete it and load directly from
-    # a zarr store for a real reconstruction run
-    data_folder = os.path.join(os.getcwd(), "data_temp")
-    os.makedirs(data_folder, exist_ok=True)
-    url = "https://zenodo.org/record/6983916/files/recOrder_test_data.zip?download=1"
-    zip_file = "recOrder_test_Data.zip"
-    output = os.path.join(data_folder, zip_file)
-    if not os.path.exists(output):
-        print("Downloading test files...")
-        download(url, out=output)
-        shutil.unpack_archive(output, extract_dir=data_folder)
-
-    temp_path = os.path.join(data_folder, timestamp + "_temp.zarr")
-    converter = TIFFConverter(
-        os.path.join(
-            data_folder,
-            "2022_08_04_recOrder_pytest_20x_04NA/2T_3P_16Z_128Y_256X_Kazansky_1",
-        ),
-        temp_path,
-    )
-    converter.run()
-
-    # setup input reader
-    reader = open_ome_zarr(temp_path, mode="r+")
-    T, C, Z, Y, X = reader["0/0/0"].data.shape
-    dtype = reader["0/0/0"].data.dtype
-    P = len(list(reader.positions()))
-    bg_data = load_bg(
-        os.path.join(data_folder, "2022_08_04_recOrder_pytest_20x_04NA/BG/"),
-        height=Y,
-        width=X,
-    )
-
-    # setup reconstructor.
-    reconstructor_args = {
-        "image_dim": (Y, X),
-        "mag": 20,  # magnification
-        "pixel_size_um": 6.5,  # pixel size in um
-        "n_slices": Z,  # number of slices in z-stack
-        "z_step_um": 2,  # z-step size in um
-        "wavelength_nm": 532,
-        "swing": 0.1,
-        "calibration_scheme": "4-State",  # "4-State" or "5-State"
-        "NA_obj": 0.4,  # numerical aperture of objective
-        "NA_illu": 0.2,  # numerical aperture of condenser
-        "n_obj_media": 1.0,  # refractive index of objective immersion media
-        "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-        "bg_correction": "global",  # BG correction method: "None", "local_fit", "global"
-        "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-        "use_gpu": False,
-        "gpu_id": 0,
-    }
-    reconstructor = initialize_reconstructor(
-        pipeline="QLIPP", **reconstructor_args
-    )
-
-    # reconstruct background Stokes
-    bg_stokes = reconstruct_qlipp_stokes(bg_data, reconstructor)
-
-    # setup output zarr
-    writer = open_ome_zarr(
-        "./output/reconstructions_" + timestamp + ".zarr",
-        layout="hcs",
-        mode="w-",
-        channel_names=[
-            "Retardance",
-            "Orientation",
-            "BF - computed",
-            "DoP",
-            "Phase",
-        ],
-    )
-    for p in range(P):
-        position = writer.create_position("0", str(p), "0")
-        position.create_zeros(name="0", shape=(T, 5, Z, Y, X), dtype=dtype)
-
-    return reader, writer, reconstructor, bg_stokes
-
-
-# define the work to be done on a single process
-# in this example a single process will loop through a set of CZYX volumes, and
-# for each volume it will read the data, apply a reconstruction, then save the
-# result to a hcs-format zarr
-def single_process(
-    reader, writer, reconstructor, bg_stokes, vol_start, vol_end
-):
-    T, C, Z, Y, X = reader["0/0/0"].data.shape
-    for vol in range(vol_start, vol_end + 1):
-        p = int(np.floor(vol / T))
-        t = int(vol % T)
-        print(f"Reconstructing vol={vol}, pos={p}, time={t}")
-
-        # read from input zarr
-        data = reader["0/" + str(p) + "/0"]["0"][t]
-
-        # reconstruct data Stokes w/ background correction
-        stokes = reconstruct_qlipp_stokes(data, reconstructor, bg_stokes)
-
-        # reconstruct birefringence from Stokes
-        birefringence = reconstruct_qlipp_birefringence(stokes, reconstructor)
-        birefringence[0] = (
-            birefringence[0] / (2 * np.pi) * (reconstructor.lambda_illu / 1000)
-        )
-
-        # reconstruct phase3D from S0
-        S0 = birefringence[2]
-        phase3D = reconstruct_phase3D(
-            S0, reconstructor, method="Tikhonov", reg_re=1e-2
-        )
-
-        # write output
-        writer["0/" + str(p) + "/0"]["0"][t, 0:4] = birefringence
-        writer["0/" + str(p) + "/0"]["0"][t, 4] = phase3D
-
-
-# prepare the processes
-# we will apply the same reconstruction to each position and time point, so we
-# need to split the CZYX volumes among the processes
-if __name__ == "__main__":
-    reader, writer, reconstructor, bg_stokes = precomputation()
-    processes = []
-    V = len(list(reader.positions())) * reader["0/0/0"].data.shape[0]
-    vol_per_process = int(np.ceil(V / N_processes))
-    print("vol_per_process", vol_per_process)
-    for i in range(N_processes):
-        vol_start = i * vol_per_process
-        vol_end = np.minimum(vol_start + vol_per_process - 1, V - 1)
-        print(f"Preparing volumes {vol_start}-{vol_end} for process {i}")
-        processes.append(
-            mp.Process(
-                target=single_process,
-                args=(
-                    reader,
-                    writer,
-                    reconstructor,
-                    bg_stokes,
-                    vol_start,
-                    vol_end,
-                ),
-            )
-        )
-    print("Starting processes...")
-    # run the processes
-    for process in processes:
-        process.start()
-    for process in processes:
-        process.join()
-
-    writer.close()
diff --git a/examples/phase.yml b/examples/phase.yml
new file mode 100644
index 00000000..d48c6aa6
--- /dev/null
+++ b/examples/phase.yml
@@ -0,0 +1,18 @@
+input_channel_names:
+- BF
+reconstruction_dimension: 3
+phase:
+  transfer_function:
+    wavelength_illumination: 0.532
+    yx_pixel_size: 0.325
+    z_pixel_size: 2.0
+    z_padding: 0
+    index_of_refraction_media: 1.3
+    numerical_aperture_detection: 1.2
+    numerical_aperture_illumination: 0.5
+    axial_flip: false
+  apply_inverse:
+    reconstruction_algorithm: Tikhonov
+    regularization_strength: 0.001
+    TV_rho_strength: 0.001
+    TV_iterations: 1
diff --git a/examples/sweep-regularization.py b/examples/sweep-regularization.py
deleted file mode 100644
index fca6570c..00000000
--- a/examples/sweep-regularization.py
+++ /dev/null
@@ -1,75 +0,0 @@
-from iohub import read_micromanager, open_ome_zarr
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_phase3D,
-)
-from recOrder.compute.phantoms import bf_3D_from_phantom
-from datetime import datetime
-import numpy as np
-import napari
-
-timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
-
-## Load a dataset
-
-# Option 1: use data simulated from a phantom and run this script as is.
-data = bf_3D_from_phantom()  # (Z, Y, X)
-
-# Option 2: load from file
-# reader = read_micromanager('/path/to/ome-tiffs/or/zarr/store/')
-# position, time, channel = 0, 0, 0
-# data = reader.get_array(position)[time, channel, ...]  # read 3D volume
-
-## Set up a reconstructor.
-Z, Y, X = data.shape
-reconstructor_args = {
-    "image_dim": (Y, X),
-    "mag": 20,  # magnification
-    "pixel_size_um": 6.5,  # pixel size in um
-    "n_slices": Z,  # number of slices in z-stack
-    "z_step_um": 2,  # z-step size in um
-    "wavelength_nm": 532,
-    "NA_obj": 0.4,  # numerical aperture of objective
-    "NA_illu": 0.2,  # numerical aperture of condenser
-    "n_obj_media": 1.0,  # refractive index of objective immersion media
-    "pad_z": 5,  # slices to pad for phase reconstruction boundary artifacts
-    "mode": "3D",  # phase reconstruction mode, "2D" or "3D"
-    "use_gpu": False,
-    "gpu_id": 0,
-}
-reconstructor = initialize_reconstructor(
-    pipeline="PhaseFromBF", **reconstructor_args
-)
-
-# Setup a single writer and viewer
-reg_powers = 10.0 ** np.arange(-3, 3)
-
-dataset = open_ome_zarr(
-    "./output/reconstructions_" + timestamp + ".zarr",
-    layout="fov",
-    mode="w-",
-    channel_names=[f"{reg:.1e}" for reg in reg_powers],
-)
-
-dataset.create_zeros(
-    "0", shape=(1, len(reg_powers), Z, Y, X), dtype=np.float32
-)
-
-v = napari.Viewer()
-v.add_image(data)
-
-# Loop through regularizations
-for i, reg in enumerate(reg_powers):
-    print(f"Reconstructing with 3D phase with reg = {reg}")
-
-    phase3D = reconstruct_phase3D(
-        data, reconstructor, method="Tikhonov", reg_re=reg
-    )
-
-    # Save each regularization into a "position" of the output zarr
-    dataset["0"][0, i, ...] = phase3D
-
-    # Add the reconstructions to the viewer
-    v.add_image(phase3D, name=f"reg = {reg}")
-
-napari.run()
diff --git a/recOrder/acq/acquisition_workers.py b/recOrder/acq/acquisition_workers.py
index 25dd073f..91bbe938 100644
--- a/recOrder/acq/acquisition_workers.py
+++ b/recOrder/acq/acquisition_workers.py
@@ -4,31 +4,25 @@
 from qtpy.QtCore import Signal
 from iohub import open_ome_zarr
 from iohub.convert import TIFFConverter
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_qlipp_birefringence,
-    reconstruct_qlipp_stokes,
-    reconstruct_phase2D,
-    reconstruct_phase3D,
+from recOrder.cli import settings
+from recOrder.cli.compute_transfer_function import (
+    compute_transfer_function_cli,
+)
+from recOrder.cli.apply_inverse_transfer_function import (
+    apply_inverse_transfer_function_cli,
 )
 from recOrder.acq.acq_functions import (
     generate_acq_settings,
     acquire_from_settings,
 )
-from recOrder.io.utils import load_bg, extract_reconstruction_parameters
-from recOrder.compute.reconstructions import QLIPPBirefringenceCompute
-from recOrder.io.metadata_reader import MetadataReader, get_last_metadata_file
-from recOrder.io.utils import ram_message, rec_bkg_to_wo_bkg
+from recOrder.io.utils import model_to_yaml
+from recOrder.io.utils import ram_message
 from napari.qt.threading import WorkerBaseSignals, WorkerBase
 from napari.utils.notifications import show_warning
 import logging
-import tifffile as tiff
 import numpy as np
 import os
-import zarr
 import shutil
-import time
-import glob
 
 # type hint/check
 from typing import TYPE_CHECKING
@@ -39,6 +33,81 @@
     from recOrder.calib.Calibration import QLIPP_Calibration
 
 
+def _generate_reconstruction_config_from_gui(
+    reconstruction_config_path,
+    mode,
+    calib_window,
+    input_channel_names,
+):
+    if mode == "birefringence" or mode == "all":
+        if calib_window.bg_option == "None":
+            background_path = ""
+            remove_estimated_background = False
+        elif calib_window.bg_option == "Measured":
+            background_path = calib_window.acq_bg_directory
+            remove_estimated_background = False
+        elif calib_window.bg_option == "Estimated":
+            background_path = ""
+            remove_estimated_background = True
+        elif calib_window.bg_option == "Measured + Estimated":
+            background_path = calib_window.acq_bg_directory
+            remove_estimated_background = True
+
+        birefringence_transfer_function_settings = (
+            settings.BirefringenceTransferFunctionSettings(
+                swing=calib_window.swing,
+            )
+        )
+        birefringence_apply_inverse_settings = settings.BirefringenceApplyInverseSettings(
+            wavelength_illumination=calib_window.recon_wavelength / 1000,
+            background_path=background_path,
+            remove_estimated_background=remove_estimated_background,
+            orientation_flip=False,  # TODO connect this to a recOrder button
+            orientation_rotate=False,  # TODO connect this to a recOrder button
+        )
+        birefringence_settings = settings.BirefringenceSettings(
+            transfer_function=birefringence_transfer_function_settings,
+            apply_inverse=birefringence_apply_inverse_settings,
+        )
+    else:
+        birefringence_settings = None
+
+    if mode == "phase" or mode == "all":
+        phase_transfer_function_settings = (
+            settings.PhaseTransferFunctionSettings(
+                wavelength_illumination=calib_window.recon_wavelength / 1000,
+                yx_pixel_size=calib_window.ps / calib_window.mag,
+                z_pixel_size=calib_window.z_step,
+                z_padding=calib_window.pad_z,
+                index_of_refraction_media=calib_window.n_media,
+                numerical_aperture_detection=calib_window.obj_na,
+                numerical_aperture_illumination=calib_window.cond_na,
+                axial_flip=False,  # TODO: connect this to a recOrder button
+            )
+        )
+        phase_apply_inverse_settings = settings.FourierApplyInverseSettings(
+            reconstruction_algorithm=calib_window.phase_regularizer,
+            regularization_strength=calib_window.ui.le_phase_strength.text(),
+            TV_rho_strength=calib_window.ui.le_rho.text(),
+            TV_iterations=calib_window.ui.le_itr.text(),
+        )
+        phase_settings = settings.PhaseSettings(
+            transfer_function=phase_transfer_function_settings,
+            apply_inverse=phase_apply_inverse_settings,
+        )
+    else:
+        phase_settings = None
+
+    reconstruction_settings = settings.ReconstructionSettings(
+        input_channel_names=input_channel_names,
+        reconstruction_dimension=int(calib_window.acq_mode[0]),
+        birefringence=birefringence_settings,
+        phase=phase_settings,
+    )
+
+    model_to_yaml(reconstruction_settings, reconstruction_config_path)
+
+
 class PolarizationAcquisitionSignals(WorkerBaseSignals):
     """
     Custom Signals class that includes napari native signals
@@ -60,16 +129,6 @@ class BFAcquisitionSignals(WorkerBaseSignals):
     aborted = Signal()
 
 
-class ListeningSignals(WorkerBaseSignals):
-    """
-    Custom Signals class that includes napari native signals
-    """
-
-    store_emitter = Signal(object)
-    dim_emitter = Signal(tuple)
-    aborted = Signal()
-
-
 class BFAcquisitionWorker(WorkerBase):
     """
     Class to execute a brightfield acquisition.  First step is to snap the images follow by a second
@@ -162,6 +221,17 @@ def work(self):
                     channel_group = group
                     break
 
+        # Create and validate reconstruction settings
+        self.config_path = os.path.join(
+            self.snap_dir, "reconstruction_settings.yml"
+        )
+        _generate_reconstruction_config_from_gui(
+            self.config_path,
+            "phase",
+            self.calib_window,
+            input_channel_names=["BF"],
+        )
+
         # Acquire 3D stack
         logging.debug("Acquiring 3D stack")
 
@@ -196,14 +266,14 @@ def work(self):
         # Reconstruct snapped images
         self.n_slices = stack.shape[2]
 
-        phase, meta = self._reconstruct(stack[0])
+        phase = self._reconstruct()
         self._check_abort()
 
-        # Save images
-        logging.debug("Saving Images")
-        self._save_imgs(phase, meta)
+        # # Save images
+        # logging.debug("Saving Images")
+        # self._save_imgs(phase, meta)
 
-        self._check_abort()
+        # self._check_abort()
 
         logging.info("Finished Acquisition")
         logging.debug("Finished Acquisition")
@@ -211,221 +281,41 @@ def work(self):
         # Emit the images and let thread know function is finished
         self.phase_image_emitter.emit(phase)
 
-    def _reconstruct(self, stack):
+    def _reconstruct(self):
         """
-        Method to reconstruct, given a 2D or 3D stack.
-        This function also checks to see if the reconstructor needs to be updated from previous acquisitions
-
-        Parameters
-        ----------
-        stack:          (nd-array) Dimensions are (C, Z, Y, X)
-
-        Returns
-        -------
-
+        Method to reconstruct
         """
-
-        self.img_dim = (stack.shape[-2], stack.shape[-1], stack.shape[-3])
-        self._check_abort()
-
-        # Initialize the reconstuctor
-
-        # if no reconstructor has been initialized before, create new reconstructor
-        if not self.calib_window.phase_reconstructor:
-            logging.debug("Computing new reconstructor")
-
-            recon = initialize_reconstructor(
-                "PhaseFromBF",
-                image_dim=(stack.shape[-2], stack.shape[-1]),
-                wavelength_nm=self.calib_window.recon_wavelength,
-                NA_obj=self.calib_window.obj_na,
-                NA_illu=self.calib_window.cond_na,
-                mag=self.calib_window.mag,
-                n_slices=self.img_dim[-1],
-                z_step_um=self.calib_window.z_step,
-                pad_z=self.calib_window.pad_z,
-                pixel_size_um=self.calib_window.ps,
-                n_obj_media=self.calib_window.n_media,
-                mode=self.calib_window.acq_mode,
-                use_gpu=self.calib_window.use_gpu,
-                gpu_id=self.calib_window.gpu_id,
-            )
-
-            # Emit reconstructor to be saved for later reconstructions
-            self.phase_reconstructor_emitter.emit(recon)
-
-        # if previous reconstructor exists
-        else:
-            self._check_abort()
-
-            # compute new reconstructor if the old reconstructor properties have been modified
-            if self._reconstructor_changed(stack.shape):
-                logging.debug(
-                    "Reconstruction settings changed, updating reconstructor"
-                )
-
-                recon = initialize_reconstructor(
-                    "PhaseFromBF",
-                    image_dim=(stack.shape[-2], stack.shape[-1]),
-                    wavelength_nm=self.calib_window.recon_wavelength,
-                    NA_obj=self.calib_window.obj_na,
-                    NA_illu=self.calib_window.cond_na,
-                    mag=self.calib_window.mag,
-                    n_slices=self.n_slices,
-                    z_step_um=self.calib_window.z_step,
-                    pad_z=self.calib_window.pad_z,
-                    pixel_size_um=self.calib_window.ps,
-                    n_obj_media=self.calib_window.n_media,
-                    mode=self.calib_window.acq_mode,
-                    use_gpu=self.calib_window.use_gpu,
-                    gpu_id=self.calib_window.gpu_id,
-                )
-
-                # Emit reconstructor to be saved for later reconstructions
-                self.phase_reconstructor_emitter.emit(recon)
-
-            # use previous reconstructor
-            else:
-                logging.debug("Using previous reconstruction settings")
-                recon = self.calib_window.phase_reconstructor
-
-        # Begin reconstruction with stokes (needed for birefringence or phase)
-        logging.debug("Reconstructing...")
         self._check_abort()
 
-        regularizer = self.calib_window.phase_regularizer
-        reg = float(self.calib_window.ui.le_phase_strength.text())
-
-        # Perform deconvolution
-        if self.dim == "2D":
-            phase = reconstruct_phase2D(
-                stack[0],
-                recon,
-                method=regularizer,
-                reg_p=reg,
-                itr=int(self.calib_window.ui.le_itr.text()),
-                rho=float(self.calib_window.ui.le_rho.text()),
-            )
-        else:
-            phase = reconstruct_phase3D(
-                stack[0],
-                recon,
-                method=regularizer,
-                reg_re=reg,
-                itr=int(self.calib_window.ui.le_itr.text()),
-                rho=float(self.calib_window.ui.le_rho.text()),
-            )
-
-        self._check_abort()
-
-        # Update metadata in zarr attributes with reconstruction parameters
-        meta = extract_reconstruction_parameters(
-            recon, magnification=self.calib_window.mag
+        # Create i/o paths
+        transfer_function_path = os.path.join(
+            self.snap_dir, "transfer_function.zarr"
+        )
+        reconstruction_path = os.path.join(
+            self.snap_dir, "reconstruction.zarr"
         )
-        meta["regularization_method"] = regularizer
-        meta["regularization_strength"] = reg
-        if regularizer == "TV":
-            meta["rho"] = float(self.calib_window.ui.le_rho.text())
-            meta["itr"] = int(self.calib_window.ui.le_itr.text())
-
-        # return both variables, could contain images or could be null
-        return phase, meta
-
-    def _save_imgs(self, phase, meta=None):
-        """
-        function to save images.
 
-        Parameters
-        ----------
-        phase:      (nd-array or None) phase image or stack
+        input_data_path = os.path.join(self.latest_out_path, "0", "0", "0")
 
-        Returns
-        -------
-
-        """
-        prefix = self.calib_window.save_name
-        name = (
-            f"ReconstructionSnap.zarr"
-            if not prefix
-            else f"{prefix}_ReconstructionSnap.zarr"
+        # TODO: skip if config files match
+        compute_transfer_function_cli(
+            input_data_path,
+            self.config_path,
+            transfer_function_path,
         )
 
-        with open_ome_zarr(
-            os.path.join(self.snap_dir, name),
-            layout="fov",
-            mode="w-",
-            channel_names=["Phase" + str(phase.ndim) + "D"],
-        ) as dataset:
-            dataset["0"] = phase[
-                (5 - phase.ndim) * (np.newaxis,) + (Ellipsis,)
-            ]
-            dataset.zattrs["recOrder"] = meta
-
-    def _reconstructor_changed(self, stack_shape: tuple):
-        """Function to check if the reconstructor has changed from the previous one in memory.
-        Serves to check if the worker attributes and reconstructor attributes have diverged.
+        apply_inverse_transfer_function_cli(
+            input_data_path,
+            transfer_function_path,
+            self.config_path,
+            reconstruction_path,
+        )
 
-        Parameters
-        ----------
-        stack_shape : tuple
-            shape of the stack array to reconstruct
+        # Read reconstruction to pass to emitters
+        with open_ome_zarr(reconstruction_path, mode="r") as dataset:
+            phase = dataset["0"][0]
 
-        Returns
-        -------
-        bool
-            whether a new reconstructor should be initialized
-        """
-        # Attributes that are directly equivalent to worker attributes
-        attr_list = {
-            "acq_mode": "phase_deconv",
-            "pad_z": "pad_z",
-            "n_media": "n_media",
-            "use_gpu": "use_gpu",
-            "gpu_id": "gpu_id",
-        }
-        # attributes that are modified upon passing them to reconstructor
-        attr_modified_list = {
-            "obj_na": "NA_obj",
-            "cond_na": "NA_illu",
-            "wavelength": "lambda_illu",
-            "n_slices": "N_defocus",
-            "swing": "chi",
-            "ps": "ps",
-        }
-        self._check_abort()
-        old_recon = self.calib_window.phase_reconstructor
-        # check if X/Y shape has changed
-        if (old_recon.N, old_recon.M) != stack_shape[-2:]:
-            return True
-        # check if equivalent attributes have diverged
-        for key, value in attr_list.items():
-            if getattr(self.calib_window, key) != getattr(old_recon, value):
-                return True
-        # modify attributes to be equivalent and check for divergence
-        for key, value in attr_modified_list.items():
-            if key == "wavelength":
-                if (
-                    self.calib_window.wavelength
-                    * 1e-3
-                    / self.calib_window.n_media
-                    != old_recon.lambda_illu
-                ):
-                    return True
-            elif key == "n_slices":
-                if getattr(self, key) != getattr(old_recon, value):
-                    return True
-            elif key == "ps":
-                if getattr(self.calib_window, key) / float(
-                    self.calib_window.mag
-                ) != getattr(old_recon, value):
-                    return True
-            else:
-                if getattr(
-                    self.calib_window, key
-                ) / self.calib_window.n_media != getattr(old_recon, value):
-                    return True
-        # not changed
-        return False
+        return phase
 
     def _cleanup_acq(self):
         # Get display windows
@@ -459,10 +349,10 @@ def _cleanup_acq(self):
                         if not save_prefix
                         else f"{save_prefix}_RawBFDataSnap.zarr"
                     )
-                    out_path = os.path.join(self.snap_dir, name)
+                    self.latest_out_path = os.path.join(self.snap_dir, name)
                     converter = TIFFConverter(
                         os.path.join(dir_, prefix),
-                        out_path,
+                        self.latest_out_path,
                         data_type="ometiff",
                         grid_layout=False,
                         label_positions=False,
@@ -565,6 +455,17 @@ def work(self):
 
         self._check_abort()
 
+        # Create and validate reconstruction settings
+        self.config_path = os.path.join(
+            self.snap_dir, "reconstruction_settings.yml"
+        )
+        _generate_reconstruction_config_from_gui(
+            self.config_path,
+            self.mode,
+            self.calib_window,
+            input_channel_names=channels,
+        )
+
         # Acquire 2D stack
         if self.dim == "2D":
             logging.debug("Acquiring 2D stack")
@@ -615,17 +516,17 @@ def work(self):
         # Reconstruct snapped images
         self._check_abort()
         self.n_slices = stack.shape[2]
-        birefringence, phase, meta = self._reconstruct(stack[0])
+        birefringence, phase = self._reconstruct()
         self._check_abort()
 
         if self.calib_window.orientation_offset:
             birefringence = self._orientation_offset(birefringence)
 
         # Save images
-        logging.debug("Saving Images")
-        self._save_imgs(birefringence, phase, meta)
+        # logging.debug("Saving Images")
+        # self._save_imgs(birefringence, phase)
 
-        self._check_abort()
+        # self._check_abort()
 
         logging.info("Finished Acquisition")
         logging.debug("Finished Acquisition")
@@ -653,7 +554,16 @@ def _orientation_offset(self, birefringence: np.ndarray):
                 "This will affect both visualization and the saved reconstruction."
             )
         )
+
+        # TODO: Move these to waveorder, connect to the `orientation_flip`
+        # and `orientation_rotate` parameters, and connect to recOrder buttons.
+
+        # 90 degree flip
         birefringence[1] = np.fmod(birefringence[1] + np.pi / 2, np.pi)
+
+        # Reflection
+        birefringence[1] = np.fmod(-birefringence[1], np.pi) + np.pi
+
         return birefringence
 
     def _check_exposure(self) -> None:
@@ -700,407 +610,49 @@ def _acquire(self) -> np.ndarray:
 
         return stack
 
-    def _reconstruct(self, stack):
+    def _reconstruct(self):
         """
-        Method to reconstruct, given a 2D or 3D stack.  First need to initialize the reconstructor given
+        Method to reconstruct.  First need to initialize the reconstructor given
         what type of acquisition it is (birefringence only skips a lot of heavy compute needed for phase).
         This function also checks to see if the reconstructor needs to be updated from previous acquisitions
 
-        Parameters
-        ----------
-        stack:          (nd-array) Dimensions are (C, Z, Y, X)
-
-        Returns
-        -------
-
         """
-
-        # get rid of z-dimension if 2D acquisition
-        stack = stack[:, 0] if self.n_slices == 1 else stack
-
         self._check_abort()
 
-        wo_background_correction = rec_bkg_to_wo_bkg(
-            self.calib_window.bg_option
+        # Create config and i/o paths
+        transfer_function_path = os.path.join(
+            self.snap_dir, "transfer_function.zarr"
         )
-
-        # Initialize the heavy reconstuctor
-        if self.mode == "phase" or self.mode == "all":
-            self._check_abort()
-
-            # if no reconstructor has been initialized before, create new reconstructor
-            if not self.calib_window.phase_reconstructor:
-                logging.debug("Computing new reconstructor")
-
-                recon = initialize_reconstructor(
-                    "QLIPP",
-                    image_dim=(stack.shape[-2], stack.shape[-1]),
-                    wavelength_nm=self.calib_window.recon_wavelength,
-                    swing=self.calib.swing,
-                    calibration_scheme=self.calib.calib_scheme,
-                    NA_obj=self.calib_window.obj_na,
-                    NA_illu=self.calib_window.cond_na,
-                    mag=self.calib_window.mag,
-                    n_slices=self.n_slices,
-                    z_step_um=self.calib_window.z_step,
-                    pad_z=self.calib_window.pad_z,
-                    pixel_size_um=self.calib_window.ps,
-                    bg_correction=wo_background_correction,
-                    n_obj_media=self.calib_window.n_media,
-                    mode=self.calib_window.acq_mode,
-                    use_gpu=self.calib_window.use_gpu,
-                    gpu_id=self.calib_window.gpu_id,
-                )
-
-                # Emit reconstructor to be saved for later reconstructions
-                self.phase_reconstructor_emitter.emit(recon)
-
-            # if previous reconstructor exists
-            else:
-                self._check_abort()
-
-                # compute new reconstructor if the old reconstructor properties have been modified
-                if self._reconstructor_changed():
-                    logging.debug(
-                        "Reconstruction settings changed, updating reconstructor"
-                    )
-
-                    recon = initialize_reconstructor(
-                        "QLIPP",
-                        image_dim=(stack.shape[-2], stack.shape[-1]),
-                        wavelength_nm=self.calib_window.recon_wavelength,
-                        swing=self.calib.swing,
-                        calibration_scheme=self.calib.calib_scheme,
-                        NA_obj=self.calib_window.obj_na,
-                        NA_illu=self.calib_window.cond_na,
-                        mag=self.calib_window.mag,
-                        n_slices=self.n_slices,
-                        z_step_um=self.calib_window.z_step,
-                        pad_z=self.calib_window.pad_z,
-                        pixel_size_um=self.calib_window.ps,
-                        bg_correction=wo_background_correction,
-                        n_obj_media=self.calib_window.n_media,
-                        mode=self.calib_window.acq_mode,
-                        use_gpu=self.calib_window.use_gpu,
-                        gpu_id=self.calib_window.gpu_id,
-                    )
-
-                    self.phase_reconstructor_emitter.emit(recon)
-
-                # use previous reconstructor
-                else:
-                    logging.debug("Using previous reconstruction settings")
-                    recon = self.calib_window.phase_reconstructor
-
-        # if phase isn't desired, initialize the lighter birefringence only reconstructor
-        # no need to save this reconstructor for later as it is pretty quick to compute
-        else:
-            self._check_abort()
-            logging.debug("Creating birefringence only reconstructor")
-            recon = initialize_reconstructor(
-                "birefringence",
-                image_dim=(stack.shape[-2], stack.shape[-1]),
-                calibration_scheme=self.calib.calib_scheme,
-                wavelength_nm=self.calib_window.recon_wavelength,
-                swing=self.calib.swing,
-                bg_correction=wo_background_correction,
-                n_slices=self.n_slices,
-            )
-
-        # Prepare background corrections for waveorder
-        # This block mimics qlipp_pipeline.py L110-119.
-        if self.calib_window.bg_option in ["global", "local_fit+"]:
-            logging.debug("Loading BG Data")
-            self._check_abort()
-            bg_data = self._load_bg(
-                self.calib_window.acq_bg_directory,
-                stack.shape[-2],
-                stack.shape[-1],
-            )
-            self._check_abort()
-            bg_stokes = recon.Stokes_recon(bg_data)
-            self._check_abort()
-            bg_stokes = recon.Stokes_transform(bg_stokes)
-            self._check_abort()
-        elif self.calib_window.bg_option == "local_fit":
-            bg_stokes = np.zeros((5, stack.shape[-2], stack.shape[-1]))
-            bg_stokes[
-                0, ...
-            ] = 1  # Set background to "identity" Stokes parameters.
-        else:
-            logging.debug("No Background Correction method chosen")
-            bg_stokes = None
-
-        # Begin reconstruction with stokes (needed for birefringence or phase)
-        logging.debug("Reconstructing...")
-        self._check_abort()
-        stokes = reconstruct_qlipp_stokes(stack, recon, bg_stokes)
-        self._check_abort()
-
-        # initialize empty variables to pass along
-        birefringence = None
-        phase = None
-
-        regularizer = self.calib_window.phase_regularizer
-
-        reg = float(self.calib_window.ui.le_phase_strength.text())
-
-        # reconstruct both phase and birefringence
-        if self.mode == "all":
-            if self.calib_window.acq_mode == "2D":
-                birefringence = reconstruct_qlipp_birefringence(
-                    np.mean(stokes, axis=1), recon
-                )
-            else:
-                birefringence = reconstruct_qlipp_birefringence(stokes, recon)
-            birefringence[0] = (
-                birefringence[0]
-                / (2 * np.pi)
-                * self.calib_window.recon_wavelength
-            )
-            self._check_abort()
-
-            if self.calib_window.acq_mode == "2D":
-                phase = reconstruct_phase2D(
-                    stokes[0],
-                    recon,
-                    method=regularizer,
-                    reg_p=reg,
-                    itr=int(self.calib_window.ui.le_itr.text()),
-                    rho=float(self.calib_window.ui.le_rho.text()),
-                )
-            else:
-                phase = reconstruct_phase3D(
-                    stokes[0],
-                    recon,
-                    method=regularizer,
-                    reg_re=reg,
-                    itr=int(self.calib_window.ui.le_itr.text()),
-                    rho=float(self.calib_window.ui.le_rho.text()),
-                )
-
-            self._check_abort()
-
-        # reconstruct phase only
-        elif self.mode == "phase":
-            if self.calib_window.acq_mode == "2D":
-                phase = reconstruct_phase2D(
-                    stokes[0],
-                    recon,
-                    method=regularizer,
-                    reg_p=reg,
-                    itr=int(self.calib_window.ui.le_itr.text()),
-                    rho=float(self.calib_window.ui.le_rho.text()),
-                )
-            else:
-                phase = reconstruct_phase3D(
-                    stokes[0],
-                    recon,
-                    method=regularizer,
-                    reg_re=reg,
-                    itr=int(self.calib_window.ui.le_itr.text()),
-                    rho=float(self.calib_window.ui.le_rho.text()),
-                )
-            self._check_abort()
-
-        # reconstruct birefringence only
-        elif self.mode == "birefringence":
-            birefringence = reconstruct_qlipp_birefringence(stokes, recon)
-            birefringence[0] = (
-                birefringence[0]
-                / (2 * np.pi)
-                * self.calib_window.recon_wavelength
-            )
-            self._check_abort()
-
-        else:
-            raise ValueError("Reconstruction Mode Not Understood")
-
-        meta = extract_reconstruction_parameters(recon, self.calib_window.mag)
-        meta["regularization_method"] = regularizer
-        meta["regularization_strength"] = reg
-        if regularizer == "TV":
-            meta["rho"] = float(self.calib_window.ui.le_rho.text())
-            meta["itr"] = int(self.calib_window.ui.le_itr.text())
-
-        # return both variables, could contain images or could be null
-        return birefringence, phase, meta
-
-    def _save_imgs(self, birefringence, phase, meta=None):
-        """
-        function to save images.
-        Makes sure file names do not overlap, i.e. nothing is overwritten.
-
-        Parameters
-        ----------
-        birefringence:      (nd-array) birefringence image(s)
-        phase:              (nd-array or None) phase image(s)
-
-        Returns
-        -------
-
-        """
-        prefix = self.calib_window.save_name
-
-        name = (
-            f"ReconstructionSnap.zarr"
-            if not prefix
-            else f"{prefix}_ReconstructionSnap.zarr"
+        reconstruction_path = os.path.join(
+            self.snap_dir, "reconstruction.zarr"
         )
 
-        with open_ome_zarr(
-            os.path.join(self.snap_dir, name),
-            layout="fov",
-            mode="w-",
-            channel_names=["Retardance", "Orientation", "BF", "Pol"],
-        ) as dataset:
-            if birefringence.ndim == 3:
-                birefringence = birefringence[:, np.newaxis]  # CYX -> CZYX
-            birefringence = birefringence[np.newaxis]  # CZYX -> TCZYX
-            dataset["0"] = birefringence
-
-            if phase is not None:
-                dataset.append_channel(
-                    "Phase" + str(phase.ndim) + "D", resize_arrays=True
-                )
-                if phase.ndim == 2:
-                    phase = phase[np.newaxis]  # YX -> ZYX
-                dataset["0"][0, 4] = phase
-
-            dataset.zattrs["recOrder"] = meta
-
-    def _load_bg(self, path, height, width):
-        """
-        # TODO: remove ROI for 1.0.0
-
-        Load background and calibration metadata.
-
-        Parameters
-        ----------
-        path:           (str) path to the BG folder
-        height:         (int) height of BG image
-        width:          (int) widht of BG image
-
-        Returns
-        -------
-
-        """
-
-        # TODO: Change to just accept ROI
-        try:
-            metadata_path = get_last_metadata_file(path)
-            metadata = MetadataReader(metadata_path)
-            roi = metadata.ROI
-        except:
-            roi = None
-
-        bg_data = load_bg(path, height, width, roi)
-
-        return bg_data
+        input_data_path = os.path.join(self.latest_out_path, "0", "0", "0")
 
-    def _reconstructor_changed(self):
-        """
-        Function to check if the reconstructor has changed from the previous one in memory.
-        Serves to check if the worker attributes and reconstructor attributes have diverged.
-
-        Returns
-        -------
+        # TODO: skip if config files match
+        compute_transfer_function_cli(
+            input_data_path,
+            self.config_path,
+            transfer_function_path,
+        )
 
-        """
+        apply_inverse_transfer_function_cli(
+            input_data_path,
+            transfer_function_path,
+            self.config_path,
+            reconstruction_path,
+        )
 
-        changed = None
-
-        # Attributes that are directly equivalent to worker attributes
-        attr_list = {
-            "acq_mode": "phase_deconv",
-            "pad_z": "pad_z",
-            "n_media": "n_media",
-            "bg_option": "bg_option",
-            "use_gpu": "use_gpu",
-            "gpu_id": "gpu_id",
-        }
-
-        # attributes that are modified upon passing them to reconstructor
-        attr_modified_list = {
-            "obj_na": "NA_obj",
-            "cond_na": "NA_illu",
-            "wavelength": "lambda_illu",
-            "n_slices": "N_defocus",
-            "swing": "chi",
-            "ps": "ps",
-        }
+        # Read reconstruction to pass to emitters
+        with open_ome_zarr(reconstruction_path, mode="r") as dataset:
+            czyx_data = dataset["0"][0]
+            birefringence = czyx_data[0:4]
+            try:
+                phase = czyx_data[4]
+            except:
+                phase = None
 
-        self._check_abort()
-        # check if equivalent attributes have diverged
-        for key, value in attr_list.items():
-            if getattr(self.calib_window, key) != getattr(
-                self.calib_window.phase_reconstructor, value
-            ):
-                changed = True
-                break
-            else:
-                changed = False
-
-        if not changed:
-            # modify attributes to be equivalent and check for divergence
-            for key, value in attr_modified_list.items():
-                if key == "swing":
-                    if self.calib.calib_scheme == "5-State":
-                        if (
-                            self.calib.swing * 2 * np.pi
-                            != self.calib_window.phase_reconstructor.chi
-                        ):
-                            changed = True
-                        else:
-                            changed = False
-                    else:
-                        if (
-                            self.calib.swing
-                            != self.calib_window.phase_reconstructor.chi
-                        ):
-                            changed = True
-                        else:
-                            changed = False
-
-                elif key == "wavelength":
-                    if (
-                        self.calib_window.recon_wavelength
-                        * 1e-3
-                        / self.calib_window.n_media
-                        != self.calib_window.phase_reconstructor.lambda_illu
-                    ):
-                        changed = True
-                        break
-                    else:
-                        changed = False
-                elif key == "n_slices":
-                    if getattr(self, key) != getattr(
-                        self.calib_window.phase_reconstructor, value
-                    ):
-                        changed = True
-                        break
-                    else:
-                        changed = False
-
-                elif key == "ps":
-                    if getattr(self.calib_window, key) / float(
-                        self.calib_window.mag
-                    ) != getattr(self.calib_window.phase_reconstructor, value):
-                        changed = True
-                        break
-                    else:
-                        changed = False
-                else:
-                    if getattr(
-                        self.calib_window, key
-                    ) / self.calib_window.n_media != getattr(
-                        self.calib_window.phase_reconstructor, value
-                    ):
-                        changed = True
-                    else:
-                        changed = False
-
-        return changed
+        return birefringence, phase
 
     def _cleanup_acq(self):
         # Get display windows
@@ -1134,10 +686,10 @@ def _cleanup_acq(self):
                         if not save_prefix
                         else f"{save_prefix}_RawPolDataSnap.zarr"
                     )
-                    out_path = os.path.join(self.snap_dir, name)
+                    self.latest_out_path = os.path.join(self.snap_dir, name)
                     converter = TIFFConverter(
                         os.path.join(dir_, prefix),
-                        out_path,
+                        self.latest_out_path,
                         data_type="ometiff",
                         grid_layout=False,
                         label_positions=False,
@@ -1149,378 +701,3 @@ def _cleanup_acq(self):
                 break
             else:
                 continue
-
-
-class ListeningWorker(WorkerBase):
-    """
-    Class to execute a birefringence/phase acquisition.  First step is to snap the images follow by a second
-    step of reconstructing those images.
-    """
-
-    def __init__(self, calib_window, bg_data):
-        super().__init__(SignalsClass=ListeningSignals)
-
-        # Save current state of GUI window
-        self.calib_window = calib_window
-
-        # Init properties
-        self.n_slices = None
-        self.n_channels = None
-        self.n_frames = None
-        self.n_pos = None
-        self.shape = None
-        self.dtype = None
-        self.root = None
-        self.prefix = None
-        self.store = None
-        self.save_directory = None
-        self.bg_data = bg_data
-        self.reconstructor = None
-
-    def _check_abort(self):
-        if self.abort_requested:
-            self.aborted.emit()
-            raise TimeoutError("Stop Requested")
-
-    def get_byte_offset(self, offsets, page):
-        """
-        Gets the byte offset from the tiff tag metadata.
-
-        210 accounts for header data + page header data.
-        162 accounts page header data
-
-        Parameters
-        ----------
-        offsets:             (dict) Offset dictionary list
-        page:               (int) Page to look at for the offset
-
-        Returns
-        -------
-        byte offset:        (int) byte offset for the image array
-
-        """
-
-        if page == 0:
-            array_offset = offsets[page] + 210
-        else:
-            array_offset = offsets[page] + 162
-
-        return array_offset
-
-    def listen_for_images(
-        self, array, file, offsets, interval, dim3, dim2, dim1, dim0, dim_order
-    ):
-        """
-
-        Parameters
-        ----------
-        array:          (nd-array) numpy array of size (C, Z)
-        file:           (string) filepath corresponding to the desired tiff image
-        offsets:        (dict) dictionary of offsets corresponding to byte offsets of pixel data in tiff image
-        interval:       (int) time interval between timepoints in seconds
-        dim3:           (int) outermost dimension to value begin at
-        dim2:           (int) first-inner dimension value to begin at
-        dim1:           (int) second-inner dimension value to begin at
-        dim0:           (int) innermost dimension value to begin at
-        dim_order:      (int) 1, 2, 3, or 4 corresponding to the dimensionality ordering of the acquisition (MM-provided)
-
-        Returns
-        -------
-        array:                      (nd-array) partially filled array of size (C, Z) to continue filling in next iteration
-        index:                      (int) current page number at end of function
-        dim3:                       (int) dimension values corresponding to where the next iteration should begin
-        dim2:                       (int) dimension values corresponding to where the next iteration should begin
-        dim1:                       (int) dimension values corresponding to where the next iteration should begin
-        dim0:                       (int) dimension values corresponding to where the next iteration should begin
-
-        """
-
-        # Order dimensions that we will loop through in order to match the acquisition
-        if dim_order == 0:
-            dims = [
-                [dim3, self.n_frames],
-                [dim2, self.n_pos],
-                [dim1, self.n_slices],
-                [dim0, self.n_channels],
-            ]
-            channel_dim = 0
-        elif dim_order == 1:
-            dims = [
-                [dim3, self.n_frames],
-                [dim2, self.n_pos],
-                [dim1, self.n_channels],
-                [dim0, self.n_slices],
-            ]
-            channel_dim = 1
-        elif dim_order == 2:
-            dims = [
-                [dim3, self.n_pos],
-                [dim2, self.n_frames],
-                [dim1, self.n_slices],
-                [dim0, self.n_channels],
-            ]
-            channel_dim = 0
-        else:
-            dims = [
-                [dim3, self.n_pos],
-                [dim2, self.n_frames],
-                [dim1, self.n_channels],
-                [dim0, self.n_slices],
-            ]
-            channel_dim = 1
-
-        # Loop through dimensions in the order they are acquired
-        idx = 0
-        for dim3 in range(dims[0][0], dims[0][1]):
-            for dim2 in range(dims[1][0], dims[1][1]):
-                for dim1 in range(dims[2][0], dims[2][1]):
-                    for dim0 in range(dims[3][0], dims[3][1]):
-                        # GET OFFSET AND WAIT
-                        if idx > 0:
-                            try:
-                                offset = self.get_byte_offset(offsets, idx)
-                            except IndexError:
-                                # NEED TO STOP BECAUSE RAN OUT OF PAGES OR REACHED END OF ACQ
-                                return array, idx, dim3, dim2, dim1, dim0
-
-                            # Checks if the offset in metadata is 0, but technically self.get_byte_offset() adds
-                            # 162 to every offset to account for tiff file header bytes
-                            while offset == 162:
-                                self._check_abort()
-                                tf = tiff.TiffFile(file)
-                                time.sleep(interval)
-                                offsets = tf.micromanager_metadata["IndexMap"][
-                                    "Offset"
-                                ]
-                                tf.close()
-                                offset = self.get_byte_offset(offsets, idx)
-                        else:
-                            offset = self.get_byte_offset(offsets, idx)
-
-                        if idx < (
-                            self.n_slices
-                            * self.n_channels
-                            * self.n_frames
-                            * self.n_pos
-                        ):
-                            # Assign dimensions based off acquisition order to correctly add image to array
-                            if dim_order == 0:
-                                t, p, c, z = dim3, dim2, dim0, dim1
-                            elif dim_order == 1:
-                                t, p, c, z = dim3, dim2, dim1, dim0
-                            elif dim_order == 2:
-                                t, p, c, z = dim2, dim3, dim0, dim1
-                            else:
-                                t, p, c, z = dim2, dim3, dim1, dim0
-
-                            self._check_abort()
-
-                            # Add Image to array
-                            img = np.memmap(
-                                file,
-                                dtype=self.dtype,
-                                mode="r",
-                                offset=offset,
-                                shape=self.shape,
-                            )
-                            array[c, z] = img
-
-                            # If Channel first, compute birefringence here
-                            if (
-                                channel_dim == 0
-                                and dim0 == self.n_channels - 1
-                            ):
-                                self._check_abort()
-
-                                # Compute birefringence
-                                self.compute_and_save(array[:, z], p, t, z)
-
-                            idx += 1
-
-                    # Reset Range to 0 to account for starting this function in middle of a dimension
-                    if idx < (
-                        self.n_slices
-                        * self.n_channels
-                        * self.n_frames
-                        * self.n_pos
-                    ):
-                        dims[2][0] = 0
-                        dims[3][0] = 0
-
-                    # If z-first, compute the birefringence here
-                    if channel_dim == 1 and dim1 == self.n_channels - 1:
-                        # Assign dimensions based off acquisition order to correctly add image to array
-                        if dim_order == 0:
-                            t, p, c, z = dim3, dim2, dim0, dim1
-                        elif dim_order == 1:
-                            t, p, c, z = dim3, dim2, dim1, dim0
-                        elif dim_order == 2:
-                            t, p, c, z = dim2, dim3, dim0, dim1
-                        else:
-                            t, p, c, z = dim2, dim3, dim1, dim0
-
-                        self._check_abort()
-                        self.compute_and_save(array, p, t, dim0)
-                    else:
-                        continue
-
-                # Reset range to 0 to account for starting this function in the middle of a dimension
-                if idx < (
-                    self.n_slices
-                    * self.n_channels
-                    * self.n_frames
-                    * self.n_pos
-                ):
-                    dims[1][0] = 0
-
-        # Return at the end of the acquisition
-        return array, idx, dim3, dim2, dim1, dim0
-
-    def compute_and_save(self, array, p, t, z):
-        if self.n_slices == 1:
-            array = array[:, 0]
-
-        birefringence = self.reconstructor.reconstruct(array)
-
-        # If the napari layer doesn't exist yet, create the zarr store to emit to napari
-        if self.prefix not in self.calib_window.viewer.layers:
-            self.store = zarr.open(
-                os.path.join(self.root, self.prefix + ".zarr")
-            )
-            self.store.zeros(
-                name="Birefringence",
-                shape=(
-                    self.n_pos,
-                    self.n_frames,
-                    2,
-                    self.n_slices,
-                    self.shape[0],
-                    self.shape[1],
-                ),
-                chunks=(1, 1, 1, 1, self.shape[0], self.shape[1]),
-                overwrite=True,
-            )
-
-            # Add data and send out the store.  Once the store has been emitted, we can add data to the store
-            # without needing to emit the store again (thanks to napari's handling of zarr datasets)
-            if len(birefringence.shape) == 4:
-                self.store["Birefringence"][p, t] = birefringence
-            else:
-                self.store["Birefringence"][p, t, :, z] = birefringence
-
-            self.store_emitter.emit(self.store)
-
-        else:
-            if len(birefringence.shape) == 4:
-                self.store["Birefringence"][p, t] = birefringence
-            else:
-                self.store["Birefringence"][p, t, :, z] = birefringence
-
-        # Emit the current dimensions so that we can update the napari dimension slider
-        self.dim_emitter.emit((p, t, z))
-
-    def work(self):
-        acq_man = self.calib_window.mm.acquisitions()
-
-        while not acq_man.isAcquisitionRunning():
-            pass
-
-        time.sleep(1)
-
-        # Get all of the dataset dimensions, settings
-        s = acq_man.getAcquisitionSettings()
-        self.root = s.root()
-        self.prefix = s.prefix()
-        self.n_frames, self.interval = (
-            (s.numFrames(), s.intervalMs() / 1000) if s.useFrames() else (1, 1)
-        )
-        self.n_channels = s.channels().size() if s.useChannels() else 1
-        self.n_slices = s.slices().size() if s.useChannels() else 1
-        self.n_pos = (
-            1
-            if not s.usePositionList()
-            else self.calib_window.mm.getPositionListManager()
-            .getPositionList()
-            .getNumberOfPositions()
-        )
-        dim_order = s.acqOrderMode()
-
-        # Get File Path corresponding to current dataset
-        path = os.path.join(self.root, self.prefix)
-        files = [fn for fn in glob.glob(path + "*") if ".zarr" not in fn]
-        index = max([int(x.strip(path + "_")) for x in files])
-
-        self.prefix = self.prefix + f"_{index}"
-        full_path = os.path.join(self.root, self.prefix)
-        first_file_path = os.path.join(
-            full_path, self.prefix + "_MMStack.ome.tif"
-        )
-
-        file = tiff.TiffFile(first_file_path)
-        file_path = first_file_path
-        self.shape = (
-            file.micromanager_metadata["Summary"]["Height"],
-            file.micromanager_metadata["Summary"]["Width"],
-        )
-        self.dtype = file.pages[0].dtype
-        offsets = file.micromanager_metadata["IndexMap"]["Offset"]
-        file.close()
-
-        self._check_abort()
-
-        # Init Reconstruction Class
-        self.reconstructor = QLIPPBirefringenceCompute(
-            self.shape,
-            self.calib_window.calib_scheme,
-            self.calib_window.wavelength,
-            self.calib_window.swing,
-            self.n_slices,
-            self.calib_window.bg_option,
-            self.bg_data,
-        )
-
-        self._check_abort()
-
-        # initialize dimensions / array for the loop
-        idx, dim3, dim2, dim1, dim0 = 0, 0, 0, 0, 0
-        array = np.zeros(
-            (self.n_channels, self.n_slices, self.shape[0], self.shape[1])
-        )
-        file_count = 0
-        total_idx = 0
-
-        # Run until the function has collected the totality of the data
-        while total_idx < (
-            self.n_slices * self.n_channels * self.n_frames * self.n_pos
-        ):
-            self._check_abort()
-
-            # this will loop through reading images in a single file as it's being written
-            # when it has successfully loaded all of the images from the file, it'll move on to the next
-            array, idx, dim3, dim2, dim1, dim0 = self.listen_for_images(
-                array,
-                file_path,
-                offsets,
-                self.interval,
-                dim3,
-                dim2,
-                dim1,
-                dim0,
-                dim_order,
-            )
-
-            total_idx += idx
-
-            # If acquisition is not finished, grab the next file and listen for images
-            if (
-                total_idx
-                != self.n_slices * self.n_channels * self.n_frames * self.n_pos
-            ):
-                time.sleep(1)
-                file_count += 1
-                file_path = os.path.join(
-                    full_path, self.prefix + f"_MMStack_{file_count}.ome.tif"
-                )
-                file = tiff.TiffFile(file_path)
-                offsets = file.micromanager_metadata["IndexMap"]["Offset"]
-                file.close()
diff --git a/recOrder/calib/Calibration.py b/recOrder/calib/Calibration.py
index d445d9d7..d8f50fbb 100644
--- a/recOrder/calib/Calibration.py
+++ b/recOrder/calib/Calibration.py
@@ -1,8 +1,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
-import matplotlib.patches as patches
-import tifffile as tiff
 import time
+from iohub import open_ome_zarr
 from recOrder.io.core_functions import *
 from recOrder.calib.Optimization import BrentOptimizer, MinScalarOptimizer
 from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
@@ -18,7 +17,7 @@
 from datetime import datetime
 from importlib_metadata import version
 
-LC_DEVICE_NAME = "MeadowlarkLcOpenSource"
+LC_DEVICE_NAME = "MeadowlarkLC"
 
 
 class QLIPP_Calibration:
@@ -107,17 +106,18 @@ def __init__(
         # Set Mode
         # TODO: make sure LC or TriggerScope are loaded in the respective modes
         allowed_modes = ["MM-Retardance", "MM-Voltage", "DAC"]
-        assert (
-            lc_control_mode in allowed_modes
-        ), f"LC control mode must be one of {allowed_modes}"
+        if lc_control_mode not in allowed_modes:
+            raise ValueError(f"LC control mode must be one of {allowed_modes}")
         self.mode = lc_control_mode
         self.LC_DAC_conversion = 4  # convert between the input range of LCs (0-20V) and the output range of the DAC (0-5V)
 
         # Initialize calibration class
         allowed_interp_methods = ["schnoor_fit", "linear"]
-        assert (
-            interp_method in allowed_interp_methods
-        ), f"LC calibration data interpolation method must be one of {allowed_interp_methods}"
+        if interp_method not in allowed_interp_methods:
+            raise ValueError(
+                "LC calibration data interpolation method must be one of "
+                f"{allowed_interp_methods}"
+            )
         dir_path = mmc.getDeviceAdapterSearchPaths().get(
             0
         )  # MM device adapter directory
@@ -294,7 +294,6 @@ def define_lc_state(self, state, lca_retardance, lcb_retardance):
             )
 
     def opt_lc(self, x, device_property, reference, normalize=False):
-
         if isinstance(x, list) or isinstance(x, tuple):
             x = x[0]
 
@@ -321,7 +320,6 @@ def opt_lc(self, x, device_property, reference, normalize=False):
             return np.abs(mean - reference)
 
     def opt_lc_cons(self, x, device_property, reference, mode):
-
         self.set_lc(x, device_property)
         swing = (self.lca_ext - x) * self.ratio
 
@@ -374,7 +372,6 @@ def opt_lc_grid(self, a_min, a_max, b_min, b_max, step):
         # coarse search
         for lca in np.arange(a_min, a_max, step):
             for lcb in np.arange(b_min, b_max, step):
-
                 self.set_lc(lca, "LCA")
                 self.set_lc(lcb, "LCB")
 
@@ -832,7 +829,6 @@ def calculate_extinction(
         )
 
     def calc_inst_matrix(self):
-
         if self.calib_scheme == "4-State":
             chi = self.swing
             inst_mat = np.array(
@@ -872,7 +868,6 @@ def calc_inst_matrix(self):
             return inst_mat
 
     def write_metadata(self, notes=None):
-
         inst_mat = self.calc_inst_matrix()
         inst_mat = np.around(inst_mat, decimals=5).tolist()
 
@@ -990,7 +985,6 @@ def _capture_state(self, state: str, n_avg: int):
         return np.mean(imgs, axis=0)
 
     def _plot_bg_images(self, imgs):
-
         img_names = (
             ["Extinction", "0", "60", "120"]
             if len(imgs) == 4
@@ -1079,39 +1073,36 @@ def capture_bg(self, n_avg, directory):
         self.mmc.setAutoShutter(False)
         self.open_shutter()
 
-        logging.debug("Capturing Background State0")
-        state0 = self._capture_state("State0", n_avg)
-        logging.debug("Saving Background State0")
-        tiff.imsave(os.path.join(directory, "State0.tif"), state0)
-
-        logging.debug("Capturing Background State1")
-        state1 = self._capture_state("State1", n_avg)
-        logging.debug("Saving Background State1")
-        tiff.imsave(os.path.join(directory, "State1.tif"), state1)
-
-        logging.debug("Capturing Background State2")
-        state2 = self._capture_state("State2", n_avg)
-        logging.debug("Saving Background State2")
-        tiff.imsave(os.path.join(directory, "State2.tif"), state2)
-
-        logging.debug("Capturing Background State3")
-        state3 = self._capture_state("State3", n_avg)
-        logging.debug("Saving Background State3")
-        tiff.imsave(os.path.join(directory, "State3.tif"), state3)
-
-        imgs = [state0, state1, state2, state3]
-
-        if self.calib_scheme == "5-State":
-            logging.debug("Capturing Background State4")
-            state4 = self._capture_state("State4", n_avg)
-            logging.debug("Saving Background State4")
-            tiff.imsave(os.path.join(directory, "State4.tif"), state4)
-            imgs.append(state4)
+        num_states = int(self.calib_scheme[0])
+
+        # Acquire background data
+        yx_list = []
+        for channel in range(num_states):
+            logging.debug(f"Capturing Background State{channel}")
+            yx_list.append(self._capture_state(f"State{channel}", n_avg))
+            logging.debug(f"Saving Background State{channel}")
+        cyx_data = np.array(yx_list)
+
+        # Save to zarr
+        with open_ome_zarr(
+            os.path.join(directory, "background.zarr"),
+            layout="hcs",
+            mode="w",
+            channel_names=[f"State{i}" for i in range(num_states)],
+        ) as dataset:
+            position = dataset.create_position("0", "0", "0")
+            position.create_zeros(
+                name="0",
+                shape=(1, num_states, 1, cyx_data.shape[1], cyx_data.shape[2]),
+                dtype=np.float32,
+                chunks=(1, 1, 1, cyx_data.shape[1], cyx_data.shape[2]),
+            )
+            position["0"][0, :, 0] = cyx_data  # save to 1C1YX array
 
         # self._plot_bg_images(np.asarray(imgs))
         self.reset_shutter()
 
-        return np.asarray(imgs)
+        return cyx_data
 
 
 class CalibrationData:
@@ -1377,7 +1368,6 @@ def interpolate_data(self, raw_data, calib_wavelengths):
             self.curve = self.spline(x_range)
 
         elif self.wavelength > 630:
-
             new_a1_y = np.interp(x_range, x_range, spline546(x_range))
             new_a2_y = np.interp(x_range, x_range, spline630(x_range))
 
@@ -1399,7 +1389,6 @@ def interpolate_data(self, raw_data, calib_wavelengths):
             self.curve = self.spline(x_range)
 
         elif 490 < self.wavelength < 546:
-
             new_a1_y = np.interp(x_range, x_range, spline490(x_range))
             new_a2_y = np.interp(x_range, x_range, spline546(x_range))
 
@@ -1416,7 +1405,6 @@ def interpolate_data(self, raw_data, calib_wavelengths):
             self.curve = self.spline(x_range)
 
         elif 546 < self.wavelength < 630:
-
             new_a1_y = np.interp(x_range, x_range, spline546(x_range))
             new_a2_y = np.interp(x_range, x_range, spline630(x_range))
 
diff --git a/recOrder/calib/calibration_workers.py b/recOrder/calib/calibration_workers.py
index 64d704ea..5ad56cca 100644
--- a/recOrder/calib/calibration_workers.py
+++ b/recOrder/calib/calibration_workers.py
@@ -4,15 +4,17 @@
 from qtpy.QtCore import Signal
 from iohub import open_ome_zarr
 from napari.qt.threading import WorkerBaseSignals, WorkerBase, thread_worker
-from recOrder.compute.reconstructions import (
-    initialize_reconstructor,
-    reconstruct_qlipp_birefringence,
-    reconstruct_qlipp_stokes,
-)
+from recOrder.cli import settings
 from recOrder.io.core_functions import set_lc_state, snap_and_average
-from recOrder.io.utils import MockEmitter
+from recOrder.io.utils import MockEmitter, model_to_yaml
 from recOrder.calib.Calibration import LC_DEVICE_NAME
 from recOrder.io.metadata_reader import MetadataReader, get_last_metadata_file
+from recOrder.cli.compute_transfer_function import (
+    compute_transfer_function_cli,
+)
+from recOrder.cli.apply_inverse_transfer_function import (
+    apply_inverse_transfer_function_cli,
+)
 import os
 import shutil
 import numpy as np
@@ -234,7 +236,6 @@ def _calibrate_4state(self):
         self._check_abort()
 
     def _calibrate_5state(self):
-
         search_radius = np.min((self.calib.swing, 0.05))
 
         self.calib.calib_scheme = "5-State"
@@ -306,7 +307,6 @@ def __init__(self, calib_window, calib):
         super().__init__(calib_window, calib)
 
     def work(self):
-
         # Make the background folder
         bg_path = os.path.join(
             self.calib_window.directory,
@@ -315,7 +315,6 @@ def work(self):
         if not os.path.exists(bg_path):
             os.mkdir(bg_path)
         else:
-
             # increment background paths
             idx = 1
             while os.path.exists(bg_path + f"_{idx}"):
@@ -334,36 +333,60 @@ def work(self):
 
         # capture and return background images
         imgs = self.calib.capture_bg(self.calib_window.n_avg, bg_path)
-        img_dim = (imgs.shape[-2], imgs.shape[-1])
-
         self.calib_window._dump_gui_state(bg_path)
 
-        # initialize reconstructor
-        recon = initialize_reconstructor(
-            "birefringence",
-            image_dim=img_dim,
-            calibration_scheme=self.calib.calib_scheme,
-            wavelength_nm=self.calib.wavelength,
-            swing=self.calib.swing,
-            bg_correction="None",
+        # build background-specific reconstruction settings
+        reconstruction_settings = settings.ReconstructionSettings(
+            input_channel_names=[
+                f"State{i}"
+                for i in range(int(self.calib_window.calib_scheme[0]))
+            ],
+            reconstruction_dimension=2,
+            birefringence=settings.BirefringenceSettings(
+                transfer_function=settings.BirefringenceTransferFunctionSettings(
+                    swing=self.calib_window.swing
+                ),
+                apply_inverse=settings.BirefringenceApplyInverseSettings(
+                    wavelength_illumination=self.calib_window.recon_wavelength
+                    / 1000,
+                    background_path="",
+                    remove_estimated_background=False,
+                    orientation_flip=False,
+                ),
+            ),
         )
 
-        self._check_abort()
-
-        # Reconstruct birefringence from BG images
-        stokes = reconstruct_qlipp_stokes(imgs, recon, None)
-
-        self._check_abort()
+        reconstruction_config_path = os.path.join(
+            bg_path, "reconstruction_settings.yml"
+        )
+        model_to_yaml(reconstruction_settings, reconstruction_config_path)
 
-        self.birefringence = reconstruct_qlipp_birefringence(stokes, recon)
+        input_data_path = os.path.join(
+            bg_path, "background.zarr", "0", "0", "0"
+        )
+        transfer_function_path = os.path.join(
+            bg_path, "transfer_function.zarr"
+        )
+        reconstruction_path = os.path.join(bg_path, "reconstruction.zarr")
 
-        self._check_abort()
+        compute_transfer_function_cli(
+            input_data_path,
+            reconstruction_config_path,
+            transfer_function_path,
+        )
 
-        # Convert retardance to nm
-        self.retardance = (
-            self.birefringence[0] / (2 * np.pi) * self.calib.wavelength
+        apply_inverse_transfer_function_cli(
+            input_data_path,
+            transfer_function_path,
+            reconstruction_config_path,
+            reconstruction_path,
         )
 
+        # Load reconstructions from file for layers
+        with open_ome_zarr(reconstruction_path, mode="r") as dataset:
+            self.retardance = dataset["0"][0, 0, 0]
+            self.birefringence = dataset["0"][0, :, 0]
+
         # Save metadata file and emit imgs
         meta_file = os.path.join(bg_path, "calibration_metadata.txt")
         self._write_meta_file(meta_file)
@@ -385,9 +408,6 @@ def work(self):
 
         self._check_abort()
 
-        self._save_bg_recon(bg_path)
-        self._check_abort()
-
         # Emit background images + background birefringence
         self.bg_image_emitter.emit(imgs)
         self.bire_image_emitter.emit((self.retardance, self.birefringence[1]))
@@ -395,45 +415,6 @@ def work(self):
         # Emit bg path
         self.bg_path_update_emitter.emit(bg_path)
 
-    def _save_bg_recon(self, bg_path: StrOrBytesPath):
-        bg_recon_path = os.path.join(bg_path, "reconstruction")
-        # create the reconstruction directory
-        if os.path.isdir(bg_recon_path):
-            shutil.rmtree(bg_recon_path)
-        elif os.path.isfile(bg_recon_path):
-            os.remove(bg_recon_path)
-        else:
-            os.mkdir(bg_recon_path)
-        # save raw reconstruction to zarr store
-        with open_ome_zarr(
-            os.path.join(bg_recon_path, "reconstruction"),
-            layout="fov",
-            mode="w-",
-            channel_names=[
-                "Retardance",
-                "Orientation",
-                "BF - computed",
-                "DoP",
-            ],
-        ) as dataset:
-            dataset["0"] = self.birefringence[np.newaxis, :, np.newaxis, ...]
-
-        # save intensity trace visualization
-        import matplotlib.pyplot as plt
-
-        plt.imsave(
-            os.path.join(bg_recon_path, "retardance.png"),
-            self.retardance,
-            cmap="gray",
-        )
-        plt.imsave(
-            os.path.join(bg_recon_path, "orientation.png"),
-            self.birefringence[1],
-            cmap="hsv",
-            vmin=0,
-            vmax=np.pi,
-        )
-
 
 @thread_worker
 def load_calibration(calib, metadata: MetadataReader):
diff --git a/recOrder/cli/apply_inverse_transfer_function.py b/recOrder/cli/apply_inverse_transfer_function.py
new file mode 100644
index 00000000..46cac81a
--- /dev/null
+++ b/recOrder/cli/apply_inverse_transfer_function.py
@@ -0,0 +1,396 @@
+import click
+import numpy as np
+import torch
+from iohub import open_ome_zarr
+from recOrder.cli.printing import echo_headline, echo_settings
+from recOrder.cli.settings import ReconstructionSettings
+from recOrder.cli.parsing import (
+    input_data_path_argument,
+    config_path_option,
+    output_dataset_option,
+)
+from recOrder.io import utils
+from waveorder.models import (
+    inplane_oriented_thick_pol3d,
+    isotropic_thin_3d,
+    phase_thick_3d,
+    isotropic_fluorescent_thick_3d,
+)
+
+
+def _check_background_consistency(background_shape, data_shape):
+    data_cyx_shape = (data_shape[1],) + data_shape[3:]
+    if background_shape != data_cyx_shape:
+        raise ValueError(
+            f"Background shape {background_shape} does not match data shape {data_cyx_shape}"
+        )
+
+
+def apply_inverse_transfer_function_cli(
+    input_data_path, transfer_function_path, config_path, output_path
+):
+    echo_headline("Starting reconstruction...")
+
+    # Load datasets
+    transfer_function_dataset = open_ome_zarr(transfer_function_path)
+    input_dataset = open_ome_zarr(input_data_path)
+
+    # Load config file
+    settings = utils.yaml_to_model(config_path, ReconstructionSettings)
+
+    # Check input channel names
+    if not set(settings.input_channel_names).issubset(
+        input_dataset.channel_names
+    ):
+        raise ValueError(
+            f"Each of the input_channel_names = {settings.input_channel_names} in {config_path} must appear in the dataset {input_data_path} which currently contains channel_names = {input_dataset.channel_names}."
+        )
+
+    # Find channel indices
+    channel_indices = []
+    for input_channel_name in settings.input_channel_names:
+        channel_indices.append(
+            input_dataset.channel_names.index(input_channel_name)
+        )
+
+    # Load dataset shape
+    t_shape = input_dataset.data.shape[0]
+
+    # Simplify important settings names
+    recon_biref = settings.birefringence is not None
+    recon_phase = settings.phase is not None
+    recon_fluo = settings.fluorescence is not None
+    recon_dim = settings.reconstruction_dimension
+
+    # Prepare output dataset
+    channel_names = []
+    if recon_biref:
+        channel_names.append("Retardance")
+        channel_names.append("Orientation")
+        channel_names.append("BF")
+        channel_names.append("Pol")
+    if recon_phase:
+        if recon_dim == 2:
+            channel_names.append("Phase2D")
+        elif recon_dim == 3:
+            channel_names.append("Phase3D")
+    if recon_fluo:
+        fluor_name = settings.input_channel_names[0]
+        if recon_dim == 2:
+            channel_names.append(fluor_name + "2D")
+        elif recon_dim == 3:
+            channel_names.append(fluor_name + "3D")
+
+    if recon_dim == 2:
+        output_z_shape = 1
+    elif recon_dim == 3:
+        output_z_shape = input_dataset.data.shape[2]
+
+    output_shape = (
+        t_shape,
+        len(channel_names),
+        output_z_shape,
+    ) + input_dataset.data.shape[3:]
+
+    # Create output dataset
+    output_dataset = open_ome_zarr(
+        output_path, layout="fov", mode="w", channel_names=channel_names
+    )
+    output_array = output_dataset.create_zeros(
+        name="0",
+        shape=output_shape,
+        dtype=np.float32,
+        chunks=(
+            1,
+            1,
+            1,
+        )
+        + input_dataset.data.shape[3:],  # chunk by YX
+    )
+
+    # Load data
+    tczyx_uint16_numpy = input_dataset.data.oindex[:, channel_indices]
+    tczyx_data = torch.tensor(
+        np.int32(tczyx_uint16_numpy), dtype=torch.float32
+    )  # convert to np.int32 (torch doesn't accept np.uint16), then convert to tensor float32
+
+    # Prepare background dataset
+    if settings.birefringence is not None:
+        biref_inverse_dict = settings.birefringence.apply_inverse.dict()
+
+        # Resolve background path into array
+        background_path = biref_inverse_dict["background_path"]
+        biref_inverse_dict.pop("background_path")
+        if background_path != "":
+            cyx_no_sample_data = utils.load_background(background_path)
+            _check_background_consistency(
+                cyx_no_sample_data.shape, input_dataset.data.shape
+            )
+        else:
+            cyx_no_sample_data = None
+
+    # Main reconstruction logic
+    # Eight different cases [2, 3] x [biref only, phase only, biref and phase, fluorescence]
+
+    # [biref only] [2 or 3]
+    if recon_biref and (not recon_phase):
+        echo_headline("Reconstructing birefringence with settings:")
+        echo_settings(settings.birefringence.apply_inverse)
+        echo_headline("Reconstructing birefringence...")
+
+        # Load transfer function
+        intensity_to_stokes_matrix = torch.tensor(
+            transfer_function_dataset["intensity_to_stokes_matrix"][0, 0, 0]
+        )
+
+        for time_index in range(t_shape):
+            # Apply
+            reconstructed_parameters = (
+                inplane_oriented_thick_pol3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index],
+                    intensity_to_stokes_matrix,
+                    cyx_no_sample_data=cyx_no_sample_data,
+                    project_stokes_to_2d=(recon_dim == 2),
+                    **biref_inverse_dict,
+                )
+            )
+            # Save
+            for param_index, parameter in enumerate(reconstructed_parameters):
+                output_array[time_index, param_index] = parameter
+
+    # [phase only]
+    if recon_phase and (not recon_biref):
+        echo_headline("Reconstructing phase with settings:")
+        echo_settings(settings.phase.apply_inverse)
+        echo_headline("Reconstructing phase...")
+
+        # check data shapes
+        if tczyx_data.shape[1] != 1:
+            raise ValueError(
+                "You have requested a phase-only reconstruction, but the input dataset has more than one channel."
+            )
+
+        # [phase only, 2]
+        if recon_dim == 2:
+            # Load transfer functions
+            absorption_transfer_function = torch.tensor(
+                transfer_function_dataset["absorption_transfer_function"][0, 0]
+            )
+            phase_transfer_function = torch.tensor(
+                transfer_function_dataset["phase_transfer_function"][0, 0]
+            )
+
+            for time_index in range(t_shape):
+                # Apply
+                (
+                    _,
+                    yx_phase,
+                ) = isotropic_thin_3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index, 0],
+                    absorption_transfer_function,
+                    phase_transfer_function,
+                    **settings.phase.apply_inverse.dict(),
+                )
+
+                # Save
+                output_array[time_index, -1, 0] = yx_phase
+
+        # [phase only, 3]
+        elif recon_dim == 3:
+            # Load transfer functions
+            real_potential_transfer_function = torch.tensor(
+                transfer_function_dataset["real_potential_transfer_function"][
+                    0, 0
+                ]
+            )
+            imaginary_potential_transfer_function = torch.tensor(
+                transfer_function_dataset[
+                    "imaginary_potential_transfer_function"
+                ][0, 0]
+            )
+
+            # Apply
+            for time_index in range(t_shape):
+                zyx_phase = phase_thick_3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index, 0],
+                    real_potential_transfer_function,
+                    imaginary_potential_transfer_function,
+                    z_padding=settings.phase.transfer_function.z_padding,
+                    z_pixel_size=settings.phase.transfer_function.z_pixel_size,
+                    wavelength_illumination=settings.phase.transfer_function.wavelength_illumination,
+                    **settings.phase.apply_inverse.dict(),
+                )
+                # Save
+                output_array[time_index, -1] = zyx_phase
+
+    # [biref and phase]
+    if recon_biref and recon_phase:
+        echo_headline("Reconstructing phase with settings:")
+        echo_settings(settings.phase.apply_inverse)
+        echo_headline("Reconstructing birefringence with settings:")
+        echo_settings(settings.birefringence.apply_inverse)
+        echo_headline("Reconstructing...")
+
+        # Load birefringence transfer function
+        intensity_to_stokes_matrix = torch.tensor(
+            transfer_function_dataset["intensity_to_stokes_matrix"][0, 0, 0]
+        )
+
+        # [biref and phase, 2]
+        if recon_dim == 2:
+            # Load phase transfer functions
+            absorption_transfer_function = torch.tensor(
+                transfer_function_dataset["absorption_transfer_function"][0, 0]
+            )
+            phase_transfer_function = torch.tensor(
+                transfer_function_dataset["phase_transfer_function"][0, 0]
+            )
+
+            for time_index in range(t_shape):
+                # Apply
+                reconstructed_parameters_2d = inplane_oriented_thick_pol3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index],
+                    intensity_to_stokes_matrix,
+                    cyx_no_sample_data=cyx_no_sample_data,
+                    project_stokes_to_2d=True,
+                    **biref_inverse_dict,
+                )
+
+                reconstructed_parameters_3d = inplane_oriented_thick_pol3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index],
+                    intensity_to_stokes_matrix,
+                    cyx_no_sample_data=cyx_no_sample_data,
+                    project_stokes_to_2d=False,
+                    **biref_inverse_dict,
+                )
+
+                brightfield_3d = reconstructed_parameters_3d[2]
+
+                (
+                    _,
+                    yx_phase,
+                ) = isotropic_thin_3d.apply_inverse_transfer_function(
+                    brightfield_3d,
+                    absorption_transfer_function,
+                    phase_transfer_function,
+                    **settings.phase.apply_inverse.dict(),
+                )
+
+                # Save
+                for param_index, parameter in enumerate(
+                    reconstructed_parameters_2d
+                ):
+                    output_array[time_index, param_index] = parameter
+                output_array[time_index, -1, 0] = yx_phase
+
+        # [biref and phase, 3]
+        elif recon_dim == 3:
+            # Load phase transfer functions
+            intensity_to_stokes_matrix = torch.tensor(
+                transfer_function_dataset["intensity_to_stokes_matrix"][
+                    0, 0, 0
+                ]
+            )
+            # Load transfer functions
+            real_potential_transfer_function = torch.tensor(
+                transfer_function_dataset["real_potential_transfer_function"][
+                    0, 0
+                ]
+            )
+            imaginary_potential_transfer_function = torch.tensor(
+                transfer_function_dataset[
+                    "imaginary_potential_transfer_function"
+                ][0, 0]
+            )
+
+            # Apply
+            for time_index in range(t_shape):
+                reconstructed_parameters_3d = inplane_oriented_thick_pol3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index],
+                    intensity_to_stokes_matrix,
+                    cyx_no_sample_data=cyx_no_sample_data,
+                    project_stokes_to_2d=False,
+                    **biref_inverse_dict,
+                )
+
+                brightfield_3d = reconstructed_parameters_3d[2]
+
+                zyx_phase = phase_thick_3d.apply_inverse_transfer_function(
+                    brightfield_3d,
+                    real_potential_transfer_function,
+                    imaginary_potential_transfer_function,
+                    z_padding=settings.phase.transfer_function.z_padding,
+                    z_pixel_size=settings.phase.transfer_function.z_pixel_size,
+                    wavelength_illumination=settings.phase.transfer_function.wavelength_illumination,
+                    **settings.phase.apply_inverse.dict(),
+                )
+                # Save
+                for param_index, parameter in enumerate(
+                    reconstructed_parameters_3d
+                ):
+                    output_array[time_index, param_index] = parameter
+                output_array[time_index, -1] = zyx_phase
+
+    # [fluo]
+    if recon_fluo:
+        echo_headline("Reconstructing fluorescence with settings:")
+        echo_settings(settings.fluorescence.apply_inverse)
+        echo_headline("Reconstructing...")
+
+        # [fluo, 2]
+        if recon_dim == 2:
+            raise NotImplementedError
+        # [fluo, 3]
+        elif recon_dim == 3:
+            # Load transfer functions
+            optical_transfer_function = torch.tensor(
+                transfer_function_dataset["optical_transfer_function"][0, 0]
+            )
+
+            # Apply
+            for time_index in range(t_shape):
+                zyx_recon = isotropic_fluorescent_thick_3d.apply_inverse_transfer_function(
+                    tczyx_data[time_index, 0],
+                    optical_transfer_function,
+                    settings.fluorescence.transfer_function.z_padding,
+                    **settings.fluorescence.apply_inverse.dict(),
+                )
+
+                # Save
+                output_array[time_index, 0] = zyx_recon
+
+    output_dataset.zattrs["settings"] = settings.dict()
+
+    echo_headline(f"Closing {output_path}\n")
+    output_dataset.close()
+    transfer_function_dataset.close()
+    input_dataset.close()
+
+    echo_headline(
+        f"Recreate this reconstruction with:\n$ recorder apply-inv-tf {input_data_path} {transfer_function_path} -c {config_path} -o {output_path}"
+    )
+
+
+@click.command()
+@click.help_option("-h", "--help")
+@input_data_path_argument()
+@click.argument(
+    "transfer_function_path",
+    type=click.Path(exists=True),
+)
+@config_path_option()
+@output_dataset_option(default="./reconstruction.zarr")
+def apply_inv_tf(
+    input_data_path, transfer_function_path, config_path, output_path
+):
+    """
+    Apply an inverse transfer function to a dataset using a configuration file.
+
+    See /examples for example configuration files.
+
+    Example usage:\n
+    $ recorder apply-inv-tf input.zarr/0/0/0 transfer-function.zarr -c /examples/birefringence.yml -o output.zarr
+    """
+    apply_inverse_transfer_function_cli(
+        input_data_path, transfer_function_path, config_path, output_path
+    )
diff --git a/recOrder/cli/compute_transfer_function.py b/recOrder/cli/compute_transfer_function.py
new file mode 100644
index 00000000..80b740aa
--- /dev/null
+++ b/recOrder/cli/compute_transfer_function.py
@@ -0,0 +1,219 @@
+import click
+import numpy as np
+from iohub import open_ome_zarr
+from recOrder.cli.printing import echo_settings, echo_headline
+from recOrder.cli.settings import ReconstructionSettings
+from recOrder.cli.parsing import (
+    input_data_path_argument,
+    config_path_option,
+    output_dataset_option,
+)
+from recOrder.io import utils
+from waveorder.models import (
+    inplane_oriented_thick_pol3d,
+    isotropic_thin_3d,
+    phase_thick_3d,
+    isotropic_fluorescent_thick_3d,
+)
+
+
+def generate_and_save_birefringence_transfer_function(settings, dataset):
+    """Generates and saves the birefringence transfer function to the dataset, based on the settings.
+
+    Parameters
+    ----------
+    settings: ReconstructionSettings
+    dataset: NGFF Node
+        The dataset that will be updated.
+    """
+    echo_headline("Generating birefringence transfer function with settings:")
+    echo_settings(settings.birefringence.transfer_function)
+
+    # Calculate transfer functions
+    intensity_to_stokes_matrix = (
+        inplane_oriented_thick_pol3d.calculate_transfer_function(
+            scheme=str(len(settings.input_channel_names)) + "-State",
+            **settings.birefringence.transfer_function.dict(),
+        )
+    )
+    # Save
+    dataset[
+        "intensity_to_stokes_matrix"
+    ] = intensity_to_stokes_matrix.cpu().numpy()[None, None, None, ...]
+
+
+def generate_and_save_phase_transfer_function(settings, dataset, zyx_shape):
+    """Generates and saves the phase transfer function to the dataset, based on the settings.
+
+    Parameters
+    ----------
+    settings: ReconstructionSettings
+    dataset: NGFF Node
+        The dataset that will be updated.
+    zyx_shape : tuple
+        A tuple of integers specifying the input data's shape in (Z, Y, X) order
+    """
+    echo_headline("Generating phase transfer function with settings:")
+    echo_settings(settings.phase.transfer_function)
+
+    if settings.reconstruction_dimension == 2:
+        # Convert zyx_shape and z_pixel_size into yx_shape and z_position_list
+        settings_dict = settings.phase.transfer_function.dict()
+        settings_dict["yx_shape"] = [zyx_shape[1], zyx_shape[2]]
+        settings_dict["z_position_list"] = list(
+            -(np.arange(zyx_shape[0]) - zyx_shape[0] // 2)
+            * settings_dict["z_pixel_size"]
+        )
+
+        # Remove unused parameters
+        settings_dict.pop("z_pixel_size")
+        settings_dict.pop("z_padding")
+
+        # Calculate transfer functions
+        (
+            absorption_transfer_function,
+            phase_transfer_function,
+        ) = isotropic_thin_3d.calculate_transfer_function(
+            **settings_dict,
+        )
+
+        # Save
+        dataset[
+            "absorption_transfer_function"
+        ] = absorption_transfer_function.cpu().numpy()[None, None, ...]
+        dataset[
+            "phase_transfer_function"
+        ] = phase_transfer_function.cpu().numpy()[None, None, ...]
+
+    elif settings.reconstruction_dimension == 3:
+        # Calculate transfer functions
+        (
+            real_potential_transfer_function,
+            imaginary_potential_transfer_function,
+        ) = phase_thick_3d.calculate_transfer_function(
+            zyx_shape=zyx_shape,
+            **settings.phase.transfer_function.dict(),
+        )
+        # Save
+        dataset[
+            "real_potential_transfer_function"
+        ] = real_potential_transfer_function.cpu().numpy()[None, None, ...]
+        dataset[
+            "imaginary_potential_transfer_function"
+        ] = imaginary_potential_transfer_function.cpu().numpy()[
+            None, None, ...
+        ]
+
+
+def generate_and_save_fluorescence_transfer_function(
+    settings, dataset, zyx_shape
+):
+    """Generates and saves the fluorescence transfer function to the dataset, based on the settings.
+
+    Parameters
+    ----------
+    settings: ReconstructionSettings
+    dataset: NGFF Node
+        The dataset that will be updated.
+    zyx_shape : tuple
+        A tuple of integers specifying the input data's shape in (Z, Y, X) order
+    """
+    echo_headline("Generating fluorescence transfer function with settings:")
+    echo_settings(settings.fluorescence.transfer_function)
+
+    if settings.reconstruction_dimension == 2:
+        raise NotImplementedError
+    elif settings.reconstruction_dimension == 3:
+        # Calculate transfer functions
+        optical_transfer_function = (
+            isotropic_fluorescent_thick_3d.calculate_transfer_function(
+                zyx_shape=zyx_shape,
+                **settings.fluorescence.transfer_function.dict(),
+            )
+        )
+        # Save
+        dataset[
+            "optical_transfer_function"
+        ] = optical_transfer_function.cpu().numpy()[None, None, ...]
+
+
+def compute_transfer_function_cli(input_data_path, config_path, output_path):
+    """CLI command to compute the transfer function given a configuration file path
+    and a desired output path.
+
+    Parameters
+    ----------
+    input_data_path : string
+        Path to the input data file.
+    config_path : string
+        Path of the configuration file.
+    output_path : string
+        Path of the output file.
+    """
+
+    # Load config file
+    settings = utils.yaml_to_model(config_path, ReconstructionSettings)
+
+    echo_headline(
+        f"Generating transfer functions and storing in {output_path}\n"
+    )
+
+    # Read shape from input dataset
+    input_dataset = open_ome_zarr(input_data_path, layout="fov", mode="r")
+    zyx_shape = input_dataset.data.shape[
+        2:
+    ]  # only loads a single position "0"
+
+    # Check input channel names
+    if not set(settings.input_channel_names).issubset(
+        input_dataset.channel_names
+    ):
+        raise ValueError(
+            f"Each of the input_channel_names = {settings.input_channel_names} in {config_path} must appear in the dataset {input_data_path} which currently contains channel_names = {input_dataset.channel_names}."
+        )
+
+    # Prepare output dataset
+    output_dataset = open_ome_zarr(
+        output_path, layout="fov", mode="w", channel_names=["None"]
+    )
+
+    # Pass settings to appropriate calculate_transfer_function and save
+    if settings.birefringence is not None:
+        generate_and_save_birefringence_transfer_function(
+            settings, output_dataset
+        )
+    if settings.phase is not None:
+        generate_and_save_phase_transfer_function(
+            settings, output_dataset, zyx_shape
+        )
+    if settings.fluorescence is not None:
+        generate_and_save_fluorescence_transfer_function(
+            settings, output_dataset, zyx_shape
+        )
+
+    # Write settings to metadata
+    output_dataset.zattrs["settings"] = settings.dict()
+
+    echo_headline(f"Closing {output_path}\n")
+    output_dataset.close()
+
+    echo_headline(
+        f"Recreate this transfer function with:\n$ recorder compute-tf {input_data_path} -c {config_path} -o {output_path}"
+    )
+
+
+@click.command()
+@click.help_option("-h", "--help")
+@input_data_path_argument()
+@config_path_option()
+@output_dataset_option(default="./transfer-function.zarr")
+def compute_tf(input_data_path, config_path, output_path):
+    """
+    Compute a transfer function using a dataset and configuration file.
+
+    See /examples/ for example configuration files.
+
+    Example usage:\n
+    $ recorder compute-tf input.zarr/0/0/0 -c /examples/birefringence.yml -o transfer_function.zarr
+    """
+    compute_transfer_function_cli(input_data_path, config_path, output_path)
diff --git a/recOrder/cli/main.py b/recOrder/cli/main.py
new file mode 100644
index 00000000..ad39a932
--- /dev/null
+++ b/recOrder/cli/main.py
@@ -0,0 +1,16 @@
+import click
+from recOrder.cli.view import view
+from recOrder.cli.reconstruct import reconstruct
+from recOrder.cli.compute_transfer_function import compute_tf
+from recOrder.cli.apply_inverse_transfer_function import apply_inv_tf
+
+
+@click.group()
+def cli():
+    """\033[92mrecOrder: Computational Toolkit for Label-Free Imaging\033[0m\n"""
+
+
+cli.add_command(view)
+cli.add_command(reconstruct)
+cli.add_command(compute_tf)
+cli.add_command(apply_inv_tf)
diff --git a/recOrder/cli/parsing.py b/recOrder/cli/parsing.py
new file mode 100644
index 00000000..298536c0
--- /dev/null
+++ b/recOrder/cli/parsing.py
@@ -0,0 +1,54 @@
+import click
+from typing import Callable
+from iohub.ngff import open_ome_zarr, Plate
+
+
+def _validate_fov_path(
+    ctx: click.Context, opt: click.Option, value: str
+) -> None:
+    dataset = open_ome_zarr(value)
+    if isinstance(dataset, Plate):
+        raise ValueError(
+            "Please supply a single position instead of an HCS plate. Likely fix: replace 'input.zarr' with 'input.zarr/0/0/0'"
+        )
+    return value
+
+
+def input_data_path_argument() -> Callable:
+    def decorator(f: Callable) -> Callable:
+        return click.argument(
+            "input-data-path",
+            type=click.Path(exists=True),
+            callback=_validate_fov_path,
+            nargs=1,
+        )(f)
+
+    return decorator
+
+
+def config_path_option() -> Callable:
+    def decorator(f: Callable) -> Callable:
+        return click.option(
+            "--config-path", "-c", required=True, help="Path to config.yml"
+        )(f)
+
+    return decorator
+
+
+def output_dataset_option(default) -> Callable:
+    click_options = [
+        click.option(
+            "--output-path",
+            "-o",
+            default=default,
+            help="Path to output.zarr",
+        )
+    ]
+    # good place to add chunking, overwrite flag, etc
+
+    def decorator(f: Callable) -> Callable:
+        for opt in click_options:
+            f = opt(f)
+        return f
+
+    return decorator
diff --git a/recOrder/cli/printing.py b/recOrder/cli/printing.py
new file mode 100644
index 00000000..8d0f6e7a
--- /dev/null
+++ b/recOrder/cli/printing.py
@@ -0,0 +1,10 @@
+import click
+import yaml
+
+
+def echo_settings(settings):
+    click.echo(yaml.dump(settings.dict(), default_flow_style=False, sort_keys=False))
+
+
+def echo_headline(headline):
+    click.echo(click.style(headline, fg="green"))
diff --git a/recOrder/cli/reconstruct.py b/recOrder/cli/reconstruct.py
new file mode 100644
index 00000000..8599a675
--- /dev/null
+++ b/recOrder/cli/reconstruct.py
@@ -0,0 +1,44 @@
+import click
+import os
+from recOrder.cli.compute_transfer_function import (
+    compute_transfer_function_cli,
+)
+from recOrder.cli.apply_inverse_transfer_function import (
+    apply_inverse_transfer_function_cli,
+)
+from recOrder.cli.parsing import (
+    input_data_path_argument,
+    config_path_option,
+    output_dataset_option,
+)
+
+
+@click.command()
+@click.help_option("-h", "--help")
+@input_data_path_argument()
+@config_path_option()
+@output_dataset_option(default="./reconstruction.zarr")
+def reconstruct(input_data_path, config_path, output_path):
+    """Reconstruct a dataset using a configuration file. This is a
+    convenience function for a `compute-tf` call followed by a `apply-inv-tf`
+    call.
+
+    See /examples for example configuration files.
+
+    Example usage:\n
+    $ recorder reconstruct input.zarr/0/0/0 -c /examples/birefringence.yml -o output.zarr
+    """
+
+    # Handle transfer function path
+    output_directory = os.path.dirname(output_path)
+    transfer_function_path = os.path.join(
+        output_directory, "transfer_function.zarr"
+    )
+
+    # Compute transfer function and apply inverse
+    compute_transfer_function_cli(
+        input_data_path, config_path, transfer_function_path
+    )
+    apply_inverse_transfer_function_cli(
+        input_data_path, transfer_function_path, config_path, output_path
+    )
diff --git a/recOrder/cli/settings.py b/recOrder/cli/settings.py
new file mode 100644
index 00000000..051dae7a
--- /dev/null
+++ b/recOrder/cli/settings.py
@@ -0,0 +1,178 @@
+import os
+from pydantic import (
+    BaseModel,
+    Extra,
+    NonNegativeInt,
+    NonNegativeFloat,
+    PositiveInt,
+    PositiveFloat,
+    root_validator,
+    validator,
+)
+from typing import Literal, List, Optional
+
+# This file defines the configuration settings for the CLI.
+
+# Example settings files in `/examples/settings/` are autmatically generated
+# by the tests in `/tests/cli_tests/test_settings.py` - `test_generate_example_settings`.
+
+# To keep the example settings up to date, run `pytest` locally when this file changes.
+
+
+# All settings classes inherit from MyBaseModel, which forbids extra parameters to guard against typos
+class MyBaseModel(BaseModel, extra=Extra.forbid):
+    pass
+
+
+# Bottom level settings
+class WavelengthIllumination(MyBaseModel):
+    wavelength_illumination: PositiveFloat = 0.532
+
+
+class BirefringenceTransferFunctionSettings(MyBaseModel):
+    swing: float = 0.1
+
+    @validator("swing")
+    def swing_range(cls, v):
+        if v <= 0 or v >= 1.0:
+            raise ValueError(f"swing = {v} should be between 0 and 1.")
+        return v
+
+
+class BirefringenceApplyInverseSettings(WavelengthIllumination):
+    background_path: str = ""
+    remove_estimated_background: bool = False
+    orientation_flip: bool = False
+    orientation_rotate: bool = False
+
+    @validator("background_path")
+    def check_background_path(cls, v):
+        if v == "":
+            return v
+
+        raw_dir = r"{}".format(v)
+        if not os.path.isdir(raw_dir):
+            raise ValueError(f"{v} is not a existing directory")
+        return raw_dir
+
+
+class FourierTransferFunctionSettings(MyBaseModel):
+    yx_pixel_size: PositiveFloat = 6.5 / 20
+    z_pixel_size: PositiveFloat = 2.0
+    z_padding: NonNegativeInt = 0
+    index_of_refraction_media: PositiveFloat = 1.3
+    numerical_aperture_detection: PositiveFloat = 1.2
+
+    @validator("numerical_aperture_detection")
+    def na_det(cls, v, values):
+        n = values["index_of_refraction_media"]
+        if v > n:
+            raise ValueError(
+                f"numerical_aperture_detection = {v} must be less than or equal to index_of_refraction_media = {n}"
+            )
+        return v
+
+    @validator("z_pixel_size")
+    def warn_unit_consistency(cls, v, values):
+        yx_pixel_size = values["yx_pixel_size"]
+        ratio = yx_pixel_size / v
+        if ratio < 1.0 / 20 or ratio > 20:
+            raise Warning(
+                f"yx_pixel_size ({yx_pixel_size}) / z_pixel_size ({v}) = {ratio}. Did you use consistent units?"
+            )
+        return v
+
+
+class FourierApplyInverseSettings(MyBaseModel):
+    reconstruction_algorithm: Literal["Tikhonov", "TV"] = "Tikhonov"
+    regularization_strength: NonNegativeFloat = 1e-3
+    TV_rho_strength: PositiveFloat = 1e-3
+    TV_iterations: NonNegativeInt = 1
+
+
+class PhaseTransferFunctionSettings(
+    FourierTransferFunctionSettings,
+    WavelengthIllumination,
+):
+    numerical_aperture_illumination: NonNegativeFloat = 0.5
+    axial_flip: bool = False
+
+    @validator("numerical_aperture_illumination")
+    def na_ill(cls, v, values):
+        n = values.get("index_of_refraction_media")
+        if v > n:
+            raise ValueError(
+                f"numerical_aperture_illumination = {v} must be less than or equal to index_of_refraction_media = {n}"
+            )
+        return v
+
+
+class FluorescenceTransferFunctionSettings(FourierTransferFunctionSettings):
+    wavelength_emission: PositiveFloat = 0.507
+
+    @validator("wavelength_emission")
+    def warn_unit_consistency(cls, v, values):
+        yx_pixel_size = values.get("yx_pixel_size")
+        ratio = yx_pixel_size / v
+        if ratio < 1.0 / 20 or ratio > 20:
+            raise Warning(
+                f"yx_pixel_size ({yx_pixel_size}) / wavelength_illumination ({v}) = {ratio}. Did you use consistent units?"
+            )
+        return v
+
+
+# Second level settings
+class BirefringenceSettings(MyBaseModel):
+    transfer_function: BirefringenceTransferFunctionSettings = (
+        BirefringenceTransferFunctionSettings()
+    )
+    apply_inverse: BirefringenceApplyInverseSettings = (
+        BirefringenceApplyInverseSettings()
+    )
+
+
+class PhaseSettings(MyBaseModel):
+    transfer_function: PhaseTransferFunctionSettings = (
+        PhaseTransferFunctionSettings()
+    )
+    apply_inverse: FourierApplyInverseSettings = FourierApplyInverseSettings()
+
+
+class FluorescenceSettings(MyBaseModel):
+    transfer_function: FluorescenceTransferFunctionSettings = (
+        FluorescenceTransferFunctionSettings()
+    )
+    apply_inverse: FourierApplyInverseSettings = FourierApplyInverseSettings()
+
+
+# Top level settings
+class ReconstructionSettings(MyBaseModel):
+    input_channel_names: List[str] = [f"State{i}" for i in range(4)]
+    reconstruction_dimension: Literal[2, 3] = 3
+    birefringence: Optional[BirefringenceSettings]
+    phase: Optional[PhaseSettings]
+    fluorescence: Optional[FluorescenceSettings]
+
+    @root_validator(pre=False)
+    def validate_reconstruction_types(cls, values):
+        if (values.get("birefringence") or values.get("phase")) and values.get(
+            "fluorescence"
+        ) is not None:
+            raise ValueError(
+                '"fluorescence" cannot be present alongside "birefringence" or "phase". Please use one configuration file for a "fluorescence" reconstruction and another configuration file for a "birefringence" and/or "phase" reconstructions.'
+            )
+        num_channel_names = len(values.get("input_channel_names"))
+        if values.get("birefringence") is None:
+            if (
+                values.get("phase") is None
+                and values.get("fluorescence") is None
+            ):
+                raise ValueError(
+                    "Provide settings for either birefringence, phase, birefringence + phase, or fluorescence."
+                )
+            if num_channel_names != 1:
+                raise ValueError(
+                    f"{num_channel_names} channels names provided. Please provide a single channel for fluorescence/phase reconstructions."
+                )
+
+        return values
diff --git a/recOrder/scripts/cli.py b/recOrder/cli/view.py
similarity index 97%
rename from recOrder/scripts/cli.py
rename to recOrder/cli/view.py
index 5e54fe53..0461aa9a 100644
--- a/recOrder/scripts/cli.py
+++ b/recOrder/cli/view.py
@@ -1,7 +1,7 @@
 import click
 import napari
 
-v = napari.Viewer()  # open viewer right away to use on hpc
+#v = napari.Viewer()  # open viewer right away to use on hpc
 import numpy as np
 from recOrder.io.utils import ret_ori_overlay
 from iohub.reader import print_info, _infer_format
@@ -170,14 +170,7 @@ def text_overlay():
     napari.run()
 
 
-@click.group()
-def cli():
-    print(
-        "\033[92mrecOrder: Computational Toolkit for Label-Free Imaging\033[0m\n"
-    )
-
-
-@cli.command()
+@click.command()
 @click.help_option("-h", "--help")
 @click.argument("filename")
 @click.option(
diff --git a/recOrder/compute/__init__.py b/recOrder/compute/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/recOrder/compute/phantoms.py b/recOrder/compute/phantoms.py
deleted file mode 100644
index e31dcccb..00000000
--- a/recOrder/compute/phantoms.py
+++ /dev/null
@@ -1,157 +0,0 @@
-import numpy as np
-import waveorder as wo
-
-
-def pol_3D_from_phantom():
-    # Adapted from waveorder/examples/2D_QLIPP_simulation/2D_QLIPP_forward.py
-
-    N = 512  # number of pixel in y dimension
-    M = 512  # number of pixel in x dimension
-    mag = 20  # magnification
-    ps = 6.5 / mag  # effective pixel size
-    lambda_illu = 0.532  # wavelength
-    n_media = 1  # refractive index in the media
-    NA_obj = 0.4  # objective NA
-    NA_illu = 0.2  # illumination NA (condenser)
-    N_defocus = 11
-    z_defocus = (
-        np.r_[:N_defocus] - N_defocus // 2
-    ) * 2  # a set of defocus planes
-    chi = 0.1 * 2 * np.pi  # swing of Polscope analyzer
-
-    # Generate phantom
-    star, theta, _ = wo.genStarTarget(N, M)
-    phase_value = 1  # average phase in radians (optical path length)
-    phi_s = star * (phase_value + 0.15)  # slower OPL across target
-    phi_f = star * (phase_value - 0.15)  # faster OPL across target
-    mu_s = np.zeros((N, M))  # absorption
-    mu_f = mu_s.copy()
-    t_eigen = np.zeros((2, N, M), complex)  # complex specimen transmission
-    t_eigen[0] = np.exp(-mu_s + 1j * phi_s)
-    t_eigen[1] = np.exp(-mu_f + 1j * phi_f)
-    _, _, fxx, fyy = wo.gen_coordinate((N, M), ps)
-    sa = theta % np.pi  # slow axes.
-
-    # Generate source
-    Source_cont = wo.gen_Pupil(fxx, fyy, NA_illu, lambda_illu)
-    Source_discrete = wo.Source_subsample(
-        Source_cont, lambda_illu * fxx, lambda_illu * fyy, subsampled_NA=0.1
-    )
-
-    # Generate simulator
-    simulator = wo.waveorder_microscopy_simulator(
-        (N, M),
-        lambda_illu,
-        ps,
-        NA_obj,
-        NA_illu,
-        z_defocus,
-        chi,
-        n_media=n_media,
-        illu_mode="Arbitrary",
-        Source=Source_discrete,
-    )
-
-    # Simulate measurements
-    I_meas, _ = simulator.simulate_waveorder_measurements(
-        t_eigen, sa, multiprocess=False
-    )
-
-    # Add background
-    photon_count = 300000
-    ext_ratio = 10000
-    const_bg = photon_count / (0.5 * (1 - np.cos(chi))) / ext_ratio
-    I_meas_noise = (
-        np.random.poisson(I_meas / np.max(I_meas) * photon_count + const_bg)
-    ).astype("float64")
-
-    # Return data and background
-    data = I_meas_noise.transpose(0, 3, 1, 2)
-    bkg = np.ones((5, N, M)) * const_bg
-    return data, bkg
-
-
-def bf_3D_from_phantom():
-    # Adapted from waveorder/3D_PODT_phase_simulation/3D_PODT_Phase_forward.py
-
-    N = 512  # number of pixel in y dimension
-    M = 512  # number of pixel in x dimension
-    mag = 20  # magnification
-    ps = 6.5 / mag  # effective pixel size
-    lambda_illu = 0.532  # wavelength
-    n_media = 1  # refractive index in the media
-    NA_obj = 0.4  # objective NA
-    NA_illu = 0.2  # illumination NA (condenser)
-    N_defocus = 31
-    psz = 2
-    z_defocus = (
-        np.r_[:N_defocus] - N_defocus // 2
-    ) * psz  # a set of defocus planes
-    chi = 0.1  # not used
-
-    # Generate sample
-    phantom = np.zeros((N, M, N_defocus), dtype=np.complex64)
-    star, _, _ = wo.genStarTarget(N, M)
-    n_sample = 1.50
-    t_obj = np.exp(1j * 2 * np.pi * psz * (star * n_sample - n_media))
-    phantom[:, :, N_defocus // 2] = t_obj
-    phantom += n_media
-
-    # Generate source
-    _, _, fxx, fyy = wo.gen_coordinate((N, M), ps)
-    Source_cont = wo.gen_Pupil(fxx, fyy, NA_illu, lambda_illu)
-    Source_discrete = wo.Source_subsample(
-        Source_cont, lambda_illu * fxx, lambda_illu * fyy, subsampled_NA=0.1
-    )
-
-    # Simulate
-    simulator = wo.waveorder_microscopy_simulator(
-        (N, M),
-        lambda_illu,
-        ps,
-        NA_obj,
-        NA_illu,
-        z_defocus,
-        chi,
-        n_media=n_media,
-        illu_mode="Arbitrary",
-        Source=Source_discrete,
-    )
-    data = simulator.simulate_3D_scalar_measurements(phantom)
-    return data.transpose((2, 0, 1))
-
-
-def fluorescence_from_phantom():
-    # Adapted from waveorder/examples/fluorescence_deconvolution/fluorescence_deconv.ipynb
-
-    N = 512  # number of pixel in y dimension
-    M = 512  # number of pixel in x dimension
-    mag = 20  # magnification
-    ps = 6.5 / mag  # effective pixel size
-    psz = 2
-    lambda_emiss = [0.532]  # wavelength
-    n_media = 1  # refractive index in the media
-    NA_obj = 0.4  # objective NA
-    N_defocus = 11
-
-    fluor_setup = wo.fluorescence_microscopy(
-        (N, M, N_defocus),
-        lambda_emiss,
-        ps,
-        psz,
-        NA_obj,
-        n_media=n_media,
-        deconv_mode="3D-WF",
-        pad_z=10,
-        use_gpu=False,
-        gpu_id=0,
-    )
-
-    star, _, _ = wo.genStarTarget(N, M)
-
-    f = np.zeros((N, M, fluor_setup.N_defocus_3D))
-    f[:, :, fluor_setup.N_defocus_3D // 2] = star
-    F = np.fft.fftn(f)
-    H = fluor_setup.OTF_WF_3D[0, ...]
-    g = np.abs(np.fft.ifftn(F * H))
-    return g.transpose((2, 0, 1))
diff --git a/recOrder/compute/reconstructions.py b/recOrder/compute/reconstructions.py
deleted file mode 100644
index 06167ddf..00000000
--- a/recOrder/compute/reconstructions.py
+++ /dev/null
@@ -1,492 +0,0 @@
-from waveorder.waveorder_reconstructor import (
-    waveorder_microscopy,
-    fluorescence_microscopy,
-)
-import numpy as np
-import time
-
-
-def initialize_reconstructor(
-    pipeline,
-    image_dim=None,
-    wavelength_nm=None,
-    swing=None,
-    calibration_scheme=None,
-    NA_obj=1.3,
-    NA_illu=0.5,
-    mag=None,
-    n_slices=None,
-    z_step_um=None,
-    in_focus_slice=None,
-    pad_z=0,
-    pixel_size_um=None,
-    bg_correction="None",
-    n_obj_media=1.3,
-    mode="3D",
-    use_gpu=False,
-    gpu_id=0,
-):
-    """
-    Initialize the QLIPP reconstructor for downstream tasks. See tags next to parameters
-    for which parameters are needed for each pipeline
-
-    Parameters
-    ----------
-
-        pipeline          : string
-                            'birefringence', 'QLIPP', 'PhaseFromBF', 'fluorescence'
-
-        image_dim         : tuple
-                            (height, width) of images in pixels
-
-        wavelength_nm      : int
-                            wavelength of illumination in nm
-
-        swing             : float
-                            swing used for calibration in waves
-
-        calibration_scheme: str
-                            '4-State' or '5-State'
-
-        NA_obj            : float
-                            numerical aperture of the detection objective
-
-        NA_illu           : float
-                            numerical aperture of the illumination condenser
-
-        mag               : float
-                            magnification used for imaging (e.g. 20 for 20x)
-
-        n_slices          : int
-                            number of slices in the z-stack
-
-        z_step_um          : float
-                            z step size of the image space
-
-        in_focus_slice     : int or None
-                            index of the in-focus z slice used for 2D-from-3D reconstructions
-                            if None, the middle z slice is used as the in focus slice
-
-        pad_z             : float
-                            how many padding slices to add for phase computation
-
-        pixel_size_um     : float
-                            pixel size of the camera in um
-
-        bg_correction      : str
-                            'local' for estimating background with scipy uniform filter
-                            'local_fit' for estimating background with polynomial fit
-                            'None' for no BG correction
-                            'Global' for normal background subtraction with the provided background
-
-        n_obj_media        : float
-                            refractive index of the objective immersion media
-
-        mode               : str
-                            '2D' or '3D' (phase reconstruction only)
-
-        use_gpu           : bool
-                            option to use gpu or not
-
-        gpu_id            : int
-                            number refering to which gpu will be used
-
-
-        Returns
-        -------
-            reconstructor     : object
-                                reconstruction object initialized with reconstruction parameters
-
-    """
-
-    anisotropy_only = False
-
-    if NA_obj > n_obj_media:
-        raise ValueError("Parameters are non-physical: NA_obj > n_obj_media")
-    if NA_illu > n_obj_media:
-        raise ValueError("Parameters are non-physical: NA_illu > n_obj_media")
-
-    if pipeline == "QLIPP" or pipeline == "PhaseFromBF":
-
-        if not NA_obj:
-            raise ValueError("Please specify NA_obj in function parameters")
-        if not NA_illu:
-            raise ValueError("Please specify NA_illu in function parameters")
-        if not mag:
-            raise ValueError(
-                "Please specify mag (magnification) in function parameters"
-            )
-        if not wavelength_nm:
-            raise ValueError(
-                "Please specify the wavelength for reconstruction"
-            )
-        if not n_slices:
-            raise ValueError("Please specify n_slices in function parameters")
-        if not z_step_um:
-            raise ValueError("Please specify z_step_um in function parameters")
-        if not pixel_size_um:
-            raise ValueError("Please specify NA_obj in function parameters")
-        if not n_obj_media:
-            raise ValueError("Please specify NA_obj in function parameters")
-        if not image_dim:
-            raise ValueError("Please specify image_dim in function parameters")
-
-        if pipeline == "QLIPP":
-            if not calibration_scheme:
-                raise ValueError(
-                    "Please specify qlipp_scheme (calibration scheme) for QLIPP reconstruction"
-                )
-            if not swing:
-                raise ValueError("Please specify swing in function parameters")
-
-    elif pipeline == "birefringence":
-
-        anisotropy_only = True
-
-        if not calibration_scheme:
-            raise ValueError(
-                "Please specify qlipp_scheme (calibration scheme) for QLIPP reconstruction"
-            )
-        if not wavelength_nm:
-            raise ValueError(
-                "Please specify the wavelength for QLIPP reconstruction"
-            )
-        if not swing:
-            raise ValueError("Please specify swing in function parameters")
-    elif pipeline == "fluorescence":
-        if not NA_obj:
-            raise ValueError(
-                "Please specify NA_obj for fluorescence reconstruction."
-            )
-    else:
-        raise ValueError(f"Pipeline {pipeline} not understood")
-
-    # Modify user inputs to fit waveorder input requirements
-    lambda_illu = wavelength_nm / 1000 if wavelength_nm else None
-    n_defocus = n_slices if n_slices else 0
-    z_step_um = 0 if not z_step_um else z_step_um
-
-    if in_focus_slice is None:
-        in_focus_slice = n_defocus // 2
-
-    z_defocus = (
-        -(np.r_[:n_defocus] - in_focus_slice) * z_step_um
-    )  # assumes stack starts from the bottom
-    ps = pixel_size_um / mag if pixel_size_um else None
-    cali = True
-    NA_obj = 0 if not NA_obj else NA_obj
-    NA_illu = 0 if not NA_illu else NA_illu
-
-    if calibration_scheme == "4-State":
-        inst_mat = np.array(
-            [
-                [1, 0, 0, -1],
-                [1, np.sin(2 * np.pi * swing), 0, -np.cos(2 * np.pi * swing)],
-                [
-                    1,
-                    -0.5 * np.sin(2 * np.pi * swing),
-                    np.sqrt(3) * np.cos(np.pi * swing) * np.sin(np.pi * swing),
-                    -np.cos(2 * np.pi * swing),
-                ],
-                [
-                    1,
-                    -0.5 * np.sin(2 * np.pi * swing),
-                    -np.sqrt(3) / 2 * np.sin(2 * np.pi * swing),
-                    -np.cos(2 * np.pi * swing),
-                ],
-            ]
-        )
-        n_channel = 4
-
-    elif calibration_scheme == "5-State":
-        swing = swing * 2 * np.pi
-        inst_mat = None
-        n_channel = 5
-
-    elif calibration_scheme == "PhaseFromBF":
-        inst_mat = None
-        n_channel = 1
-        swing = 0
-    else:
-        inst_mat = None
-        n_channel = 1
-        swing = 0
-
-    print("Initializing Reconstructor...")
-    start_time = time.time()
-    if pipeline != "fluorescence":
-        recon = waveorder_microscopy(
-            img_dim=image_dim,
-            lambda_illu=lambda_illu,
-            ps=ps,
-            NA_obj=NA_obj,
-            NA_illu=NA_illu,
-            z_defocus=z_defocus,
-            chi=swing,
-            n_media=n_obj_media,
-            cali=cali,
-            bg_option=bg_correction,
-            A_matrix=inst_mat,
-            QLIPP_birefringence_only=anisotropy_only,
-            pad_z=pad_z,
-            phase_deconv=mode,
-            illu_mode="BF",
-            use_gpu=use_gpu,
-            gpu_id=gpu_id,
-        )
-        recon.N_channel = n_channel
-    else:
-
-        recon = fluorescence_microscopy(
-            img_dim=image_dim + (n_slices,),
-            lambda_emiss=[lambda_illu],
-            ps=ps,
-            psz=z_step_um,
-            NA_obj=NA_obj,
-            n_media=n_obj_media,
-            deconv_mode="3D-WF",
-            pad_z=pad_z,
-            use_gpu=use_gpu,
-            gpu_id=gpu_id,
-        )
-
-    elapsed_time = (time.time() - start_time) / 60
-    print(f"Finished Initializing Reconstructor ({elapsed_time:0.2f} min)")
-
-    return recon
-
-
-def reconstruct_qlipp_stokes(data, recon, bg_stokes=None):
-    """
-    From intensity data, use the waveorder.waveorder_microscopy (recon) to build a stokes array
-        if recon background correction flag is selected, will also perform backgroudn correction
-
-    Parameters
-    ----------
-        data                : np.ndarray or zarr array
-                              intensity data of shape: (C, Z, Y, X)
-
-        recon               : waveorder.waveorder_microscopy object
-                              initialized by initialize_reconstructor
-
-        bg_stokes           : np.ndarray (5, Y, X) or (4, Y, X)
-                              stokes array representing background data
-
-    Returns
-    -------
-        stokes_stack        : np.ndarray
-                              array representing stokes array
-                              has shape: (C, Z, Y, X), where C = 5
-                              or (C, Y, X) if not reconstructing a z-stack
-    """
-
-    stokes_data = recon.Stokes_recon(np.copy(data))
-    stokes_data = recon.Stokes_transform(stokes_data)
-
-    # Don't do background correction if BG data isn't provided
-    if recon.bg_option == "None" or bg_stokes is None:
-        return stokes_data  # C(Z)YX
-
-    # Compute Stokes with background correction
-    else:
-        if len(np.shape(stokes_data)) == 4:
-            s_image = recon.Polscope_bg_correction(
-                np.transpose(stokes_data, (-4, -2, -1, -3)), bg_stokes
-            )
-            s_image = np.transpose(s_image, (0, 3, 1, 2))  # Tranpose to CZYX
-        else:
-            s_image = recon.Polscope_bg_correction(
-                stokes_data, bg_stokes
-            )  # CYX
-
-        return s_image
-
-
-def reconstruct_qlipp_birefringence(stokes, recon):
-    """
-    From stokes data, use waveorder.waveorder_microscopy (recon) to build a birefringence array
-
-    Parameters
-    ----------
-    stokes                  : np.ndarray or zarr array
-                              stokes array generated by reconstruct_qlipp_stokes
-                              dimensions: (C, Z, Y, X) or (C, Y, X)
-
-    recon                   : waveorder.waveorder_microscopy object
-                              initialized by initialize_reconstructor
-
-    Returns
-    -------
-        recon_data          : np.ndarray
-                              volume of shape (C, Z, Y, X) or (C, Y, X) containing reconstructed birefringence data.
-    """
-
-    if stokes.ndim != 4 and stokes.ndim != 3:
-        raise ValueError(f"Incompatible stokes dimension: {stokes.shape}")
-
-    return recon.Polarization_recon(stokes)
-
-
-def reconstruct_phase2D(
-    S0, recon, method="Tikhonov", reg_p=1e-4, rho=1, lambda_p=1e-4, itr=50
-):
-    """
-    Reconstruct 2D phase from a given S0 or BF stack.
-
-    Parameters
-    ----------
-    S0:             (nd-array) BF/S0 stack of dimensions (Z, Y, X)
-    recon:          (waveorder_microscopy Object): initialized reconstructor object
-    method:         (str) Regularization method 'Tikhonov' or 'TV'
-    reg_p:          (float) Tikhonov regularization parameters
-    rho:            (float) TV regularization parameter
-    lambda_p:       (float) TV regularization parameter
-    itr:            (int) TV Regularization number of iterations
-
-    Returns
-    -------
-    phase2D:        (nd-array) Phase2D image of size (Y, X)
-
-    """
-
-    S0 = np.transpose(S0, (1, 2, 0))
-    _, phase2D = recon.Phase_recon(
-        np.copy(S0).astype("float"),
-        method=method,
-        reg_p=reg_p,
-        rho=rho,
-        lambda_p=lambda_p,
-        itr=itr,
-        verbose=False,
-    )
-
-    return phase2D
-
-
-def reconstruct_phase3D(
-    S0, recon, method="Tikhonov", reg_re=1e-4, rho=1e-3, lambda_re=1e-4, itr=50
-):
-    """
-    Reconstruct 2D phase from a given S0 or BF stack.
-
-    Parameters
-    ----------
-    S0:             (nd-array) BF/S0 stack of dimensions (Z, Y, X)
-    recon:          (waveorder_microscopy Object): initialized reconstructor object
-    method:         (str) Regularization method 'Tikhonov' or 'TV'
-    reg_p:          (float) Tikhonov regularization parameters
-    rho:            (float) TV regularization parameter
-    lambda_p:       (float) TV regularization parameter
-    itr:            (int) TV Regularization number of iterations
-
-    Returns
-    -------
-    phase3D:        (nd-array) Phase2D image of size (Z, Y, X)
-
-    """
-
-    S0 = np.transpose(S0, (1, 2, 0))
-    phase3D = recon.Phase_recon_3D(
-        np.copy(S0).astype("float"),
-        method=method,
-        reg_re=reg_re,
-        rho=rho,
-        lambda_re=lambda_re,
-        itr=itr,
-        verbose=False,
-    )
-
-    phase3D = np.transpose(phase3D, (-1, -3, -2))
-
-    return phase3D
-
-
-def reconstruct_density_from_fluorescence(data3D, recon, bg_level=0, reg=1e-2):
-    """
-    Reconstruct 3D density from fluorescence intensity
-
-    Parameters
-    ----------
-    data3D:         (nd-array) Stack of dimensions (Z, Y, X)
-    recon:          (fluorescence_microscopy Object): initialized reconstructor object
-    reg:            (float) Tikhonov regularization parameters
-
-    Returns
-    -------
-    density:        (nd-array) density image of size (Z, Y, X)
-
-    """
-
-    data3D = np.transpose(data3D, (1, 2, 0))
-    density = recon.deconvolve_fluor_3D(
-        np.copy(data3D).astype("float"),
-        bg_level=[bg_level],
-        reg=[reg],
-    )
-
-    return np.transpose(density, (-1, -3, -2))
-
-
-class QLIPPBirefringenceCompute:
-    """
-    Convenience Class for computing QLIPP birefringence only.
-    """
-
-    def __init__(
-        self,
-        shape,
-        scheme,
-        wavelength,
-        swing,
-        n_slices,
-        bg_option,
-        bg_data=None,
-    ):
-
-        self.shape = shape
-        self.scheme = scheme
-        self.wavelength = wavelength
-        self.swing = swing
-        self.n_slices = n_slices
-        self.bg_option = bg_option
-
-        self.reconstructor = initialize_reconstructor(
-            pipeline="birefringence",
-            image_dim=self.shape,
-            calibration_scheme=self.scheme,
-            wavelength_nm=self.wavelength,
-            swing=self.swing,
-            bg_correction=self.bg_option,
-            n_slices=self.n_slices,
-        )
-
-        if bg_option != "None":
-            self.bg_stokes = reconstruct_qlipp_stokes(
-                bg_data, self.reconstructor
-            )
-        else:
-            self.bg_stokes = None
-
-    def reconstruct(self, array):
-        """
-        reconstructs raw data into birefringence data
-
-        Parameters
-        ----------
-        array:          (nd-array) of image shape (C, Z, Y, X) or (C, Y, X) with minimum C=4
-
-        Returns
-        -------
-        birefringence:  (nd-array) birefringence array of size (2, Z, Y, X) or (2, Y, X)
-                                    first channel is retardance [nm] second channel is orientation [0, pi]
-
-        """
-        stokes = reconstruct_qlipp_stokes(
-            array, self.reconstructor, self.bg_stokes
-        )
-        birefringence = reconstruct_qlipp_birefringence(
-            stokes, self.reconstructor
-        )
-        birefringence[0] = birefringence[0] / (2 * np.pi) * self.wavelength
-
-        return birefringence[0:2]
diff --git a/recOrder/io/utils.py b/recOrder/io/utils.py
index 3c6731ab..e1b418b3 100644
--- a/recOrder/io/utils.py
+++ b/recOrder/io/utils.py
@@ -1,15 +1,21 @@
+from typing import Literal
 import glob
 import logging
 import os
 import psutil
+import torch
 import textwrap
 import tifffile as tiff
 import numpy as np
+import yaml
 from colorspacious import cspace_convert
+from iohub import open_ome_zarr
 from matplotlib.colors import hsv_to_rgb
 from waveorder.waveorder_reconstructor import waveorder_microscopy
+from recOrder.cli import settings
 
 
+# TO BE DEPRECATED
 def extract_reconstruction_parameters(reconstructor, magnification=None):
     """
     Function that extracts the reconstruction parameters from a waveorder reconstructor.  Works for waveorder_microscopy class.
@@ -53,91 +59,12 @@ def extract_reconstruction_parameters(reconstructor, magnification=None):
     return attr_dict
 
 
-def load_bg(bg_path, height, width, ROI=None):
-    """
-    Parameters
-    ----------
-    bg_path         : (str) path to the folder containing background images
-    height          : (int) height of image in pixels # Remove for 1.0.0
-    width           : (int) width of image in pixels # Remove for 1.0.0
-    ROI             : (tuple)  ROI of the background images to use, if None, full FOV will be used # Remove for 1.0.0
-
-    Returns
-    -------
-    bg_data   : (ndarray) Array of background data w/ dimensions (N_channel, Y, X)
-    """
-
-    bg_paths = glob.glob(os.path.join(bg_path, "*.tif"))
-    bg_paths.sort()
-
-    # Backwards compatibility warning
-    if ROI is not None and ROI != (
-        0,
-        0,
-        width,
-        height,
-    ):  # TODO: Remove for 1.0.0
-        warning_msg = """
-        Earlier versions of recOrder (0.1.2 and earlier) would have averaged over the background ROI. 
-        This behavior is now considered a bug, and future versions of recOrder (0.2.0 and later) 
-        will not average over the background. 
-        """
-        logging.warning(warning_msg)
-
-    # Load background images
-    bg_img_list = []
-    for bg_path in bg_paths:
-        bg_img_list.append(tiff.imread(bg_path))
-    bg_img_arr = np.array(bg_img_list)  # CYX
-
-    # Error if shapes do not match
-    # TODO: 1.0.0 move these validation check to waveorder's Polscope_bg_correction
-    if bg_img_arr.shape[1:] != (height, width):
-        error_msg = "The background image has a different X/Y size than the acquired image."
-        raise ValueError(error_msg)
-
-    return bg_img_arr  # CYX
-
-
-def create_grid_from_coordinates(xy_coords, rows, columns):
-    """
-    Function to create a grid from XY-position coordinates.  Useful for generating HCS Zarr metadata.
-
-    Parameters
-    ----------
-    xy_coords:          (list) XY Stage position list in the order in which it was acquired: (X, Y) tuple.
-    rows:               (int) number of rows in the grid-like acquisition
-    columns:            (int) number of columns in the grid-like acquisition
-
-    Returns
-    -------
-    pos_index_grid      (array) A grid-like array mimicking the shape of the acquisition where the value in the array
-                                corresponds to the position index at that location.
-    """
-
-    coords = dict()
-    coords_list = []
-    for idx, pos in enumerate(xy_coords):
-        coords[idx] = pos
-        coords_list.append(pos)
-
-    # sort by X and then by Y
-    coords_list.sort(key=lambda x: x[0])
-    coords_list.sort(key=lambda x: x[1])
-
-    # reshape XY coordinates into their proper 2D shape
-    grid = np.reshape(coords_list, (rows, columns, 2))
-    pos_index_grid = np.zeros((rows, columns), "uint16")
-    keys = list(coords.keys())
-    vals = list(coords.values())
-
-    for row in range(rows):
-        for col in range(columns):
-
-            # append position index (key) into a final grid by indexed into the coordinate map (values)
-            pos_index_grid[row, col] = keys[vals.index(list(grid[row, col]))]
-
-    return pos_index_grid
+def load_background(background_path):
+    with open_ome_zarr(
+        os.path.join(background_path, "background.zarr", "0", "0", "0")
+    ) as dataset:
+        cyx_data = dataset["0"][0, :, 0]
+        return torch.tensor(cyx_data, dtype=torch.float32)
 
 
 class MockEmitter:
@@ -145,48 +72,6 @@ def emit(self, value):
         pass
 
 
-def get_unimodal_threshold(input_image):
-    """Determines optimal unimodal threshold
-    https://users.cs.cf.ac.uk/Paul.Rosin/resources/papers/unimodal2.pdf
-    https://www.mathworks.com/matlabcentral/fileexchange/45443-rosin-thresholding
-    :param np.array input_image: generate mask for this image
-    :return float best_threshold: optimal lower threshold for the foreground
-     hist
-    """
-
-    hist_counts, bin_edges = np.histogram(
-        input_image,
-        bins=256,
-        range=(input_image.min(), np.percentile(input_image, 99.5)),
-    )
-    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
-
-    # assuming that background has the max count
-    max_idx = np.argmax(hist_counts)
-    int_with_max_count = bin_centers[max_idx]
-    p1 = [int_with_max_count, hist_counts[max_idx]]
-
-    # find last non-empty bin
-    pos_counts_idx = np.where(hist_counts > 0)[0]
-    last_binedge = pos_counts_idx[-1]
-    p2 = [bin_centers[last_binedge], hist_counts[last_binedge]]
-
-    best_threshold = -np.inf
-    max_dist = -np.inf
-    for idx in range(max_idx, last_binedge, 1):
-        x0 = bin_centers[idx]
-        y0 = hist_counts[idx]
-        a = [p1[0] - p2[0], p1[1] - p2[1]]
-        b = [x0 - p2[0], y0 - p2[1]]
-        cross_ab = a[0] * b[1] - b[0] * a[1]
-        per_dist = np.linalg.norm(cross_ab) / np.linalg.norm(a)
-        if per_dist > max_dist:
-            best_threshold = x0
-            max_dist = per_dist
-    assert best_threshold > -np.inf, "Error in unimodal thresholding"
-    return best_threshold
-
-
 def ram_message():
     """
     Determine if the system's RAM capacity is sufficient for running reconstruction.
@@ -282,7 +167,9 @@ def generic_hsv_overlay(
     return overlay_final[0] if mode == "2D" else overlay_final
 
 
-def ret_ori_overlay(retardance, orientation, ret_max=10, cmap="JCh"):
+def ret_ori_overlay(
+    retardance, orientation, ret_max=10, cmap: Literal["JCh", "HSV"] = "JCh"
+):
     """
     This function will create an overlay of retardance and orientation with two different colormap options.
     HSV is the standard Hue, Saturation, Value colormap while JCh is a similar colormap but is perceptually uniform.
@@ -301,14 +188,14 @@ def ret_ori_overlay(retardance, orientation, ret_max=10, cmap="JCh"):
     """
     if retardance.shape != orientation.shape:
         raise ValueError(
-            f"Retardance and Orientation shapes do not match: {retardance.shape} vs. {orientation.shape}"
+            "Retardance and Orientation shapes do not match: "
+            f"{retardance.shape} vs. {orientation.shape}"
         )
 
     # Prepare input and output arrays
     ret_ = np.clip(retardance, 0, ret_max)  # clip and copy
-    ori_ = (
-        np.copy(orientation) * 360 / np.pi
-    )  # convert 180 degree range into 360 to match periodicity of hue.
+    # Convert 180 degree range into 360 to match periodicity of hue.
+    ori_ = orientation * 360 / np.pi
     overlay_final = np.zeros_like(retardance)
 
     # FIX ME: this binning code leads to artifacts.
@@ -350,3 +237,96 @@ def ret_ori_overlay(retardance, orientation, ret_max=10, cmap="JCh"):
         raise ValueError(f"Colormap {cmap} not understood")
 
     return overlay_final
+
+
+def model_to_yaml(model, yaml_path):
+    """
+    Save a model's dictionary representation to a YAML file.
+
+    Parameters
+    ----------
+    model : object
+        The model object to convert to YAML.
+    yaml_path : str
+        The path to the output YAML file.
+
+    Raises
+    ------
+    TypeError
+        If the `model` object does not have a `dict()` method.
+
+    Notes
+    -----
+    This function converts a model object into a dictionary representation
+    using the `dict()` method. It removes any fields with None values before
+    writing the dictionary to a YAML file.
+
+    Examples
+    --------
+    >>> from my_model import MyModel
+    >>> model = MyModel()
+    >>> model_to_yaml(model, 'model.yaml')
+
+    """
+    if not hasattr(model, "dict"):
+        raise TypeError("The 'model' object does not have a 'dict()' method.")
+
+    model_dict = model.dict()
+
+    # Remove None-valued fields
+    clean_model_dict = {
+        key: value for key, value in model_dict.items() if value is not None
+    }
+
+    with open(yaml_path, "w+") as f:
+        yaml.dump(
+            clean_model_dict, f, default_flow_style=False, sort_keys=False
+        )
+
+
+def yaml_to_model(yaml_path, model):
+    """
+    Load model settings from a YAML file and create a model instance.
+
+    Parameters
+    ----------
+    yaml_path : str
+        The path to the YAML file containing the model settings.
+    model : class
+        The model class used to create an instance with the loaded settings.
+
+    Returns
+    -------
+    object
+        An instance of the model class with the loaded settings.
+
+    Raises
+    ------
+    TypeError
+        If the provided model is not a class or does not have a callable constructor.
+    FileNotFoundError
+        If the YAML file specified by `yaml_path` does not exist.
+
+    Notes
+    -----
+    This function loads model settings from a YAML file using `yaml.safe_load()`.
+    It then creates an instance of the provided `model` class using the loaded settings.
+
+    Examples
+    --------
+    >>> from my_model import MyModel
+    >>> model = yaml_to_model('model.yaml', MyModel)
+
+    """
+    if not callable(getattr(model, "__init__", None)):
+        raise TypeError(
+            "The provided model must be a class with a callable constructor."
+        )
+
+    try:
+        with open(yaml_path, "r") as file:
+            raw_settings = yaml.safe_load(file)
+    except FileNotFoundError:
+        raise FileNotFoundError(f"The YAML file '{yaml_path}' does not exist.")
+
+    return model(**raw_settings)
diff --git a/recOrder/plugin/main_widget.py b/recOrder/plugin/main_widget.py
index 985d3b2e..599e44c5 100644
--- a/recOrder/plugin/main_widget.py
+++ b/recOrder/plugin/main_widget.py
@@ -1,44 +1,51 @@
 # TODO: remove in Python 3.11
 from __future__ import annotations
 
-from recOrder.calib.Calibration import QLIPP_Calibration, LC_DEVICE_NAME
+import json
+import logging
+import os
+import textwrap
+import time
+from os.path import dirname
+from pathlib import Path, PurePath
+
+# type hint/check
+from typing import TYPE_CHECKING
+
+import dask.array as da
+import numpy as np
+import yaml
+from dask import delayed
+from napari import Viewer
+from napari.components import LayerList
+from napari.qt.threading import create_worker
+from napari.utils.events import Event
+from napari.utils.notifications import show_info, show_warning
+from numpy.typing import NDArray
+from numpydoc.docscrape import NumpyDocString
+from packaging import version
 from pycromanager import Core, Studio, zmq_bridge
-from qtpy.QtCore import Slot, Signal, Qt
-from qtpy.QtWidgets import QWidget, QFileDialog, QSizePolicy, QSlider
-from qtpy.QtGui import QPixmap, QColor
+from qtpy.QtCore import Qt, Signal, Slot
+from qtpy.QtGui import QColor, QPixmap
+from qtpy.QtWidgets import QFileDialog, QSizePolicy, QSlider, QWidget
 from superqt import QDoubleRangeSlider, QRangeSlider
+from waveorder.waveorder_reconstructor import waveorder_microscopy
+
+from recOrder.acq.acquisition_workers import (
+    BFAcquisitionWorker,
+    PolarizationAcquisitionWorker,
+)
 from recOrder.calib import Calibration
+from recOrder.calib.Calibration import LC_DEVICE_NAME, QLIPP_Calibration
 from recOrder.calib.calibration_workers import (
-    CalibrationWorker,
     BackgroundCaptureWorker,
+    CalibrationWorker,
     load_calibration,
 )
-from recOrder.acq.acquisition_workers import (
-    PolarizationAcquisitionWorker,
-    BFAcquisitionWorker,
-)
-from recOrder.plugin import gui
 from recOrder.io.core_functions import set_lc_state, snap_and_average
 from recOrder.io.metadata_reader import MetadataReader
 from recOrder.io.utils import ret_ori_overlay
-from waveorder.waveorder_reconstructor import waveorder_microscopy
-from pathlib import Path, PurePath
-from napari import Viewer
-from napari.utils.notifications import show_warning, show_info
-from napari.qt.threading import create_worker
-from numpydoc.docscrape import NumpyDocString
-from packaging import version
-import numpy as np
-from numpy.typing import NDArray
-import os
-from os.path import dirname
-import json
-import logging
-import textwrap
-import yaml
-
-# type hint/check
-from typing import TYPE_CHECKING
+from recOrder.plugin import gui
 
 # avoid runtime import error
 if TYPE_CHECKING:
@@ -218,6 +225,12 @@ def __init__(self, napari_viewer: Viewer):
         self.ui.qbutton_acq_ret_ori_phase.clicked[bool].connect(
             self.acq_ret_ori_phase
         )
+
+        # hook to render overlay
+        # acquistion updates existing layers and moves them to the top which triggers this event
+        self.viewer.layers.events.moved.connect(self.handle_layers_updated)
+        self.viewer.layers.events.inserted.connect(self.handle_layers_updated)
+
         # Commenting for 0.3.0. Consider debugging or deleting for 1.0.0.
         # self.ui.cb_colormap.currentIndexChanged[int].connect(
         #     self.enter_colormap
@@ -303,7 +316,7 @@ def __init__(self, napari_viewer: Viewer):
         self.orientation_offset = False
         self.pad_z = 0
         self.phase_reconstructor = None
-        self.acq_bg_directory = None
+        self.acq_bg_directory = ""
         self.auto_shutter = True
         self.lca_dac = None
         self.lcb_dac = None
@@ -786,7 +799,7 @@ def _check_requirements_for_acq(self, mode):
 
         """
         # check if a QLIPP_Calibration object has been initialized
-        if not self.calib:
+        if mode != "phase" and not self.calib:
             raise RuntimeError("Please run or load calibration first.")
 
         # initialize the variable to keep track of the success of the requirement check
@@ -817,7 +830,10 @@ def _check_requirements_for_acq(self, mode):
                 raise_error = True
 
             # check background path if 'Measured' or 'Measured + Estimated' is selected
-            if self.bg_option == "local_fit+" or self.bg_option == "global":
+            if (
+                self.bg_option == "Measured"
+                or self.bg_option == "Measured + Estimated"
+            ):
                 success = self._check_line_edit("bg_path")
                 if not success:
                     raise_error = True
@@ -879,7 +895,7 @@ def connect_to_mm(self):
         -------
 
         """
-        RECOMMENDED_MM = "20220920"
+        RECOMMENDED_MM = "20230426"
         ZMQ_TARGET_VERSION = "4.2.0"
         try:
             self.mmc = Core(convert_camel_case=False)
@@ -1127,6 +1143,7 @@ def _add_or_update_image_layer(
         if name in self.viewer.layers:
             self.viewer.layers[name].data = image
             if move_to_top:
+                logging.debug(f"Moving layer {name} to the top.")
                 src_index = self.viewer.layers.index(name)
                 self.viewer.layers.move(src_index, dest_index=-1)
         else:
@@ -1146,29 +1163,62 @@ def handle_bg_bire_image_update(self, value):
             value[1], "Background Orientation", cmap="hsv"
         )
 
-    def _draw_bire_overlay(self, overlay):
-        self._add_or_update_image_layer(
-            overlay, "BirefringenceOverlay" + self.acq_mode, cmap="rgb"
+    def handle_layers_updated(self, event: Event):
+        layers: LayerList = event.source
+        latest_layer_name = layers[-1].name
+        channels = ["Retardance", "Orientation"]
+        for ch in channels:
+            if latest_layer_name.startswith(ch):
+                suffix = latest_layer_name.replace(ch, "")
+                other_name = channels[1 - channels.index(ch)] + suffix
+                overlay_name = "BirefringenceOverlay" + suffix
+                if other_name in layers and overlay_name not in layers:
+                    logging.info(
+                        "Detected updated birefringence layers: "
+                        f"'{latest_layer_name}', '{other_name}'"
+                    )
+                    self._draw_bire_overlay(
+                        [ch + suffix, other_name], overlay_name
+                    )
+        if latest_layer_name.startswith(channels[1]):
+            logging.info(
+                "Detected orientation layer in updated layer list."
+                "Setting its colormap to HSV."
+            )
+            self.viewer.layers[latest_layer_name].colormap = "hsv"
+
+    def _draw_bire_overlay(
+        self, source_names: list[str, str], overlay_name: str
+    ):
+        def _layer_data(name: str):
+            return self.viewer.layers[name].data
+
+        def _draw(overlay):
+            self._add_or_update_image_layer(overlay, overlay_name, cmap="rgb")
+
+        orientation_name, retardance_name = sorted(source_names)
+        retardance = _layer_data(retardance_name)
+        orientation = _layer_data(orientation_name)
+        worker = create_worker(
+            ret_ori_overlay,
+            retardance=retardance,
+            orientation=orientation,
+            ret_max=np.percentile(np.ravel(retardance), 99.99),
+            cmap=self.colormap,
         )
+        worker.start()
+        logging.info(f"Updating the birefringence overlay layer.")
+        if retardance.size >= 2 * 2048 * 2048:
+            show_info("Generating large overlay. This might take a moment...")
+        worker.returned.connect(_draw)
 
     @Slot(object)
     def handle_bire_image_update(self, value: NDArray):
         # generate overlay in a separate thread
-        overlay_worker = create_worker(
-            ret_ori_overlay,
-            retardance=value[0],
-            orientation=value[1],
-            ret_max=np.percentile(value[0], 99.99),
-            cmap=self.colormap,
-        )
-        overlay_worker.returned.connect(self._draw_bire_overlay)
-        overlay_worker.start()
         for i, channel in enumerate(("Retardance", "Orientation")):
-            name = channel + self.acq_mode
+            name = channel
             cmap = "gray" if channel != "Orientation" else "hsv"
             self._add_or_update_image_layer(value[i], name, cmap=cmap)
-        if self.acq_mode == "3D":
-            show_info("Generating 3D overlay. This might take a moment...")
 
     @Slot(object)
     def handle_phase_image_update(self, value):
@@ -1521,17 +1571,17 @@ def enter_bg_correction(self):
             self.ui.label_bg_path.setHidden(False)
             self.ui.le_bg_path.setHidden(False)
             self.ui.qbutton_browse_bg_path.setHidden(False)
-            self.bg_option = "global"
+            self.bg_option = "Measured"
         elif state == 2:
             self.ui.label_bg_path.setHidden(True)
             self.ui.le_bg_path.setHidden(True)
             self.ui.qbutton_browse_bg_path.setHidden(True)
-            self.bg_option = "local_fit"
+            self.bg_option = "Estimated"
         elif state == 3:
             self.ui.label_bg_path.setHidden(False)
             self.ui.le_bg_path.setHidden(False)
             self.ui.qbutton_browse_bg_path.setHidden(False)
-            self.bg_option = "local_fit+"
+            self.bg_option = "Measured + Estimated"
 
     @Slot()
     def enter_gpu_id(self):
diff --git a/setup.cfg b/setup.cfg
index 8bf83809..d79769ee 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -37,8 +37,8 @@ python_requires = >=3.9
 setup_requires = setuptools_scm
 # add your package requirements here
 install_requires =
-	waveorder==1.0.0rc0
-	pycromanager==0.24.1
+	waveorder==2.0.0rc1
+	pycromanager==0.27.2
 	click>=8.0.1
 	natsort>=7.1.1
 	colorspacious>=1.1.2
@@ -59,14 +59,15 @@ dev =
 	tox
 	pre-commit
 	black
+	hypothesis
 
 [options.package_data]
 * = *.yaml
 
 [options.entry_points]
 console_scripts =
-	recorder = recOrder.scripts.cli:cli
-	recOrder = recOrder.scripts.cli:cli
+	recorder = recOrder.cli.main:cli
+	recOrder = recOrder.cli.main:cli
 	
 napari.manifest =
 	recOrder = recOrder:napari.yaml


Test patch:
diff --git a/.github/workflows/test_and_deploy.yml b/.github/workflows/test_and_deploy.yml
index 4c037308..89481a23 100644
--- a/.github/workflows/test_and_deploy.yml
+++ b/.github/workflows/test_and_deploy.yml
@@ -11,7 +11,7 @@ on:
       - "*"
   pull_request:
     branches:
-      - main
+      - "*"
   workflow_dispatch:
 
 jobs:
diff --git a/recOrder/tests/cli_tests/test_cli.py b/recOrder/tests/cli_tests/test_cli.py
new file mode 100644
index 00000000..e23d2180
--- /dev/null
+++ b/recOrder/tests/cli_tests/test_cli.py
@@ -0,0 +1,10 @@
+from recOrder.cli.main import cli
+from click.testing import CliRunner
+
+
+def test_main():
+    runner = CliRunner()
+    result = runner.invoke(cli)
+
+    assert result.exit_code == 0
+    assert "Toolkit" in result.output
diff --git a/recOrder/tests/cli_tests/test_compute_tf.py b/recOrder/tests/cli_tests/test_compute_tf.py
new file mode 100644
index 00000000..d4d56328
--- /dev/null
+++ b/recOrder/tests/cli_tests/test_compute_tf.py
@@ -0,0 +1,121 @@
+from recOrder.cli.main import cli
+from click.testing import CliRunner
+from recOrder.cli import settings
+from recOrder.io import utils
+from recOrder.cli.compute_transfer_function import (
+    generate_and_save_phase_transfer_function,
+    generate_and_save_birefringence_transfer_function,
+    generate_and_save_fluorescence_transfer_function,
+)
+
+
+def test_compute_transfer(tmp_path, input_zarr):
+    recon_settings = settings.ReconstructionSettings(
+        input_channel_names=[f"State{i}" for i in range(4)],
+        reconstruction_dimension=3,
+        birefringence=settings.BirefringenceSettings(),
+        phase=settings.PhaseSettings(),
+    )
+    config_path = tmp_path / "test.yml"
+    utils.model_to_yaml(recon_settings, config_path)
+
+    path, _ = input_zarr
+    runner = CliRunner()
+    result = runner.invoke(
+        cli, ["compute-tf", str(path), "-c", str(config_path)]
+    )
+    assert result.exit_code == 0
+
+
+def test_compute_transfer_blank_config():
+    runner = CliRunner()
+    for option in ("-c ", "--config-path "):
+        cmd = "compute-tf " + option
+        result = runner.invoke(cli, cmd)
+        assert result.exit_code == 2
+        assert "Error" in result.output
+
+
+def test_compute_transfer_blank_output():
+    runner = CliRunner()
+    for option in ("-o ", "--output-path "):
+        cmd = "compute-tf " + option
+        result = runner.invoke(cli, cmd)
+        assert result.exit_code == 2
+        assert "Error" in result.output
+
+
+def test_compute_transfer_output_file(tmp_path, input_zarr):
+    recon_settings = settings.ReconstructionSettings(
+        input_channel_names=["BF"],
+        reconstruction_dimension=3,
+        phase=settings.PhaseSettings(),
+    )
+    config_path = tmp_path / "test.yml"
+    utils.model_to_yaml(recon_settings, config_path)
+
+    input_path, _ = input_zarr
+    runner = CliRunner()
+    for option in ("-o", "--output-path"):
+        for output_folder in ["test1", "test2/test"]:
+            output_path = tmp_path.joinpath(output_folder)
+            result = runner.invoke(
+                cli,
+                [
+                    "compute-tf",
+                    str(input_path),
+                    "-c",
+                    str(config_path),
+                    str(option),
+                    str(output_path),
+                ],
+            )
+            assert result.exit_code == 0
+            assert str(output_path) in result.output
+            assert output_path.exists()
+
+
+def test_stokes_matrix_write(birefringence_phase_recon_settings_function):
+    settings, dataset = birefringence_phase_recon_settings_function
+    generate_and_save_birefringence_transfer_function(settings, dataset)
+    assert dataset["intensity_to_stokes_matrix"]
+
+
+def test_absorption_and_phase_write(
+    birefringence_phase_recon_settings_function,
+):
+    settings, dataset = birefringence_phase_recon_settings_function
+    generate_and_save_phase_transfer_function(settings, dataset, (3, 4, 5))
+    assert dataset["real_potential_transfer_function"]
+    assert dataset["imaginary_potential_transfer_function"]
+    assert dataset["imaginary_potential_transfer_function"].shape == (
+        1,
+        1,
+        3,
+        4,
+        5,
+    )
+    assert "absorption_transfer_function" not in dataset
+    assert "phase_transfer_function" not in dataset
+
+
+def test_phase_3dim_write(birefringence_phase_recon_settings_function):
+    settings, dataset = birefringence_phase_recon_settings_function
+    settings.reconstruction_dimension = 2
+    generate_and_save_phase_transfer_function(settings, dataset, (3, 4, 5))
+    assert dataset["absorption_transfer_function"]
+    assert dataset["phase_transfer_function"]
+    assert dataset["phase_transfer_function"].shape == (1, 1, 3, 4, 5)
+    assert "real_potential_transfer_function" not in dataset
+    assert "imaginary_potential_transfer_function" not in dataset
+
+
+def test_fluorescence_write(fluorescence_recon_settings_function):
+    settings, dataset = fluorescence_recon_settings_function
+    generate_and_save_fluorescence_transfer_function(
+        settings, dataset, (3, 4, 5)
+    )
+    assert dataset["optical_transfer_function"]
+    assert dataset["optical_transfer_function"].shape == (1, 1, 3, 4, 5)
+    assert "real_potential_transfer_function" not in dataset
+    assert "imaginary_potential_transfer_function" not in dataset
diff --git a/recOrder/tests/cli_tests/test_reconstruct.py b/recOrder/tests/cli_tests/test_reconstruct.py
new file mode 100644
index 00000000..cdf4f8c8
--- /dev/null
+++ b/recOrder/tests/cli_tests/test_reconstruct.py
@@ -0,0 +1,101 @@
+import numpy as np
+from recOrder.cli.main import cli
+from recOrder.cli import settings
+from recOrder.io import utils
+from click.testing import CliRunner
+from iohub.ngff import open_ome_zarr
+
+
+def test_reconstruct(tmp_path):
+    input_path = tmp_path / "input.zarr"
+
+    # Generate input "dataset"
+    channel_names = [f"State{x}" for x in range(4)]
+    dataset = open_ome_zarr(
+        input_path,
+        layout="fov",
+        mode="w",
+        channel_names=channel_names,
+    )
+    dataset.create_zeros("0", (2, 4, 4, 5, 6), dtype=np.uint16)
+
+    # Setup options
+    birefringence_settings = settings.BirefringenceSettings(
+        transfer_function=settings.BirefringenceTransferFunctionSettings()
+    )
+    all_options = [
+        (birefringence_settings, None, 2),
+        (birefringence_settings, settings.PhaseSettings(), 2),
+        (birefringence_settings, None, 3),
+        (birefringence_settings, settings.PhaseSettings(), 3),
+    ]
+
+    for birefringence_option, phase_option, dimension_option in all_options:
+        if (birefringence_option is None) and (phase_option is None):
+            continue
+
+        # Generate recon settings
+        recon_settings = settings.ReconstructionSettings(
+            input_channel_names=channel_names,
+            reconstruction_dimension=dimension_option,
+            birefringence=birefringence_option,
+            phase=phase_option,
+        )
+        config_path = tmp_path / "test.yml"
+        utils.model_to_yaml(recon_settings, config_path)
+
+        # Run CLI
+        runner = CliRunner()
+        tf_path = input_path.with_name("tf.zarr")
+        runner.invoke(
+            cli,
+            [
+                "compute-tf",
+                str(input_path),
+                "-c",
+                str(config_path),
+                "-o",
+                str(tf_path),
+            ],
+        )
+        assert tf_path.exists()
+
+        # Apply the tf
+        result_path = input_path.with_name("result.zarr")
+
+        result_inv = runner.invoke(
+            cli,
+            [
+                "apply-inv-tf",
+                str(input_path),
+                str(tf_path),
+                "-c",
+                str(config_path),
+                "-o",
+                str(result_path),
+            ],
+        )
+        assert result_path.exists()
+        assert result_inv.exit_code == 0
+        assert "Reconstructing" in result_inv.output
+
+        # Check output
+        result_dataset = open_ome_zarr(result_path)
+        assert result_dataset["0"].shape[0] == 2
+        assert result_dataset["0"].shape[3:] == (5, 6)
+
+        # Test direct recon
+        result_inv = runner.invoke(
+            cli,
+            [
+                "reconstruct",
+                str(input_path),
+                "-c",
+                str(config_path),
+                "-o",
+                str(result_path),
+            ],
+        )
+        assert result_path.exists()
+        assert result_inv.exit_code == 0
+        assert "Reconstructing" in result_inv.output
diff --git a/recOrder/tests/cli_tests/test_settings.py b/recOrder/tests/cli_tests/test_settings.py
new file mode 100644
index 00000000..4e37bc19
--- /dev/null
+++ b/recOrder/tests/cli_tests/test_settings.py
@@ -0,0 +1,122 @@
+import pytest
+from recOrder.cli import settings
+from recOrder.io import utils
+from pydantic import ValidationError
+
+
+def test_reconstruction_settings():
+    # Test defaults
+    s = settings.ReconstructionSettings(
+        birefringence=settings.BirefringenceSettings()
+    )
+    assert len(s.input_channel_names) == 4
+    assert s.birefringence.apply_inverse.background_path == ""
+    assert s.phase == None
+    assert s.fluorescence == None
+
+    # Test logic that "fluorescence" or ("phase" and/or "birefringence")
+    s = settings.ReconstructionSettings(
+        input_channel_names=["GFP"],
+        birefringence=None,
+        phase=None,
+        fluorescence=settings.FluorescenceSettings(),
+    )
+
+    assert s.fluorescence.apply_inverse.reconstruction_algorithm == "Tikhonov"
+
+    # Not allowed to supply both phase/biref and fluorescence
+    with pytest.raises(ValidationError):
+        settings.ReconstructionSettings(
+            phase=settings.PhaseSettings(),
+            fluorescence=settings.FluorescenceSettings(),
+        )
+
+    # Test incorrect settings
+    with pytest.raises(ValidationError):
+        settings.ReconstructionSettings(input_channel_names=3)
+
+    with pytest.raises(ValidationError):
+        settings.ReconstructionSettings(reconstruction_dimension=1)
+
+    # Test typo
+    with pytest.raises(ValidationError):
+        settings.ReconstructionSettings(
+            flurescence=settings.FluorescenceSettings()
+        )
+
+
+def test_biref_tf_settings():
+    settings.BirefringenceTransferFunctionSettings(swing=0.1)
+
+    with pytest.raises(ValidationError):
+        settings.BirefringenceTransferFunctionSettings(swing=1.1)
+
+    with pytest.raises(ValidationError):
+        settings.BirefringenceTransferFunctionSettings(scheme="Test")
+
+
+def test_phase_tf_settings():
+    settings.PhaseTransferFunctionSettings(
+        index_of_refraction_media=1.0, numerical_aperture_detection=0.8
+    )
+
+    with pytest.raises(ValidationError):
+        settings.PhaseTransferFunctionSettings(
+            index_of_refraction_media=1.0, numerical_aperture_detection=1.1
+        )
+
+    # Inconsistent units
+    with pytest.raises(Warning):
+        settings.PhaseTransferFunctionSettings(
+            yx_pixel_size=650, z_pixel_size=0.3
+        )
+
+    # Extra parameter
+    with pytest.raises(ValidationError):
+        settings.PhaseTransferFunctionSettings(zyx_pixel_size=650)
+
+
+def test_fluor_tf_settings():
+    settings.FluorescenceTransferFunctionSettings(
+        wavelength_emission=0.500, yx_pixel_size=0.2
+    )
+
+    with pytest.raises(Warning):
+        settings.FluorescenceTransferFunctionSettings(
+            wavelength_emission=0.500, yx_pixel_size=2000
+        )
+
+
+def test_generate_example_settings():
+    example_path = "./examples/"
+
+    s0 = settings.ReconstructionSettings(
+        birefringence=settings.BirefringenceSettings(),
+        phase=settings.PhaseSettings(),
+    )
+    s1 = settings.ReconstructionSettings(
+        input_channel_names=["BF"],
+        phase=settings.PhaseSettings(),
+    )
+    s2 = settings.ReconstructionSettings(
+        birefringence=settings.BirefringenceSettings(),
+    )
+    s3 = settings.ReconstructionSettings(
+        input_channel_names=["GFP"],
+        fluorescence=settings.FluorescenceSettings(),
+    )
+    file_names = [
+        "birefringence-and-phase.yml",
+        "phase.yml",
+        "birefringence.yml",
+        "fluorescence.yml",
+    ]
+    settings_list = [s0, s1, s2, s3]
+
+    # Save to examples folder and test roundtrip
+    for file_name, settings_obj in zip(file_names, settings_list):
+        utils.model_to_yaml(settings_obj, example_path + file_name)
+        settings_roundtrip = utils.yaml_to_model(
+            example_path + file_name, settings.ReconstructionSettings
+        )
+        assert settings_obj.dict() == settings_roundtrip.dict()
diff --git a/recOrder/tests/conftest.py b/recOrder/tests/conftest.py
index e69de29b..a6bd4275 100644
--- a/recOrder/tests/conftest.py
+++ b/recOrder/tests/conftest.py
@@ -0,0 +1,48 @@
+import pytest
+import numpy as np
+from recOrder.cli import settings
+from iohub.ngff import open_ome_zarr
+
+
+@pytest.fixture(scope="function")
+def input_zarr(tmp_path):
+    path = tmp_path / "input.zarr"
+
+    dataset = open_ome_zarr(
+        path,
+        layout="fov",
+        mode="w",
+        channel_names=[f"State{i}" for i in range(4)] + ["BF"],
+    )
+    dataset.create_zeros("0", (2, 5, 4, 5, 6), dtype=np.uint16)
+    yield path, dataset
+
+
+@pytest.fixture(scope="function")
+def birefringence_phase_recon_settings_function(tmp_path):
+    recon_settings = settings.ReconstructionSettings(
+        birefringence=settings.BirefringenceSettings(),
+        phase=settings.PhaseSettings(),
+    )
+    dataset = open_ome_zarr(
+        tmp_path,
+        layout="fov",
+        mode="w",
+        channel_names=[f"State{i}" for i in range(4)],
+    )
+    yield recon_settings, dataset
+
+
+@pytest.fixture(scope="function")
+def fluorescence_recon_settings_function(tmp_path):
+    recon_settings = settings.ReconstructionSettings(
+        input_channel_names=["GFP"],
+        fluorescence=settings.FluorescenceSettings(),
+    )
+    dataset = open_ome_zarr(
+        tmp_path,
+        layout="fov",
+        mode="w",
+        channel_names=[f"State{i}" for i in range(4)],
+    )
+    yield recon_settings, dataset
diff --git a/recOrder/tests/util_tests/test_io.py b/recOrder/tests/util_tests/test_io.py
new file mode 100644
index 00000000..b523cd5f
--- /dev/null
+++ b/recOrder/tests/util_tests/test_io.py
@@ -0,0 +1,47 @@
+import os
+import pytest
+import yaml
+
+from recOrder.cli import settings
+from recOrder.io.utils import model_to_yaml, yaml_to_model
+
+
+@pytest.fixture
+def model():
+    # Create a sample model object
+    return settings.ReconstructionSettings(
+        birefringence=settings.BirefringenceSettings()
+    )
+
+
+@pytest.fixture
+def yaml_path(tmpdir):
+    # Create a temporary YAML file path
+    return os.path.join(tmpdir, "model.yaml")
+
+
+def test_model_to_yaml(model, yaml_path):
+    # Call the function under test
+    model_to_yaml(model, yaml_path)
+
+    # Check if the YAML file is created
+    assert os.path.exists(yaml_path)
+
+    # Load the YAML file and verify its contents
+    with open(yaml_path, "r") as f:
+        yaml_data = yaml.safe_load(f)
+
+    # Check if the YAML data is a dictionary
+    assert isinstance(yaml_data, dict)
+
+    # Check YAML data
+    assert "input_channel_names" in yaml_data
+
+
+def test_model_to_yaml_invalid_model():
+    # Create an object that does not have a 'dict()' method
+    invalid_model = "not a model object"
+
+    # Call the function and expect a TypeError
+    with pytest.raises(TypeError):
+        model_to_yaml(invalid_model, "model.yaml")
diff --git a/recOrder/tests/util_tests/test_overlays.py b/recOrder/tests/util_tests/test_overlays.py
new file mode 100644
index 00000000..aa27f966
--- /dev/null
+++ b/recOrder/tests/util_tests/test_overlays.py
@@ -0,0 +1,66 @@
+import hypothesis.extra.numpy as npst
+import hypothesis.strategies as st
+import numpy as np
+from hypothesis import given
+from numpy.typing import NDArray
+from numpy.testing import assert_equal
+
+from recOrder.io.utils import ret_ori_overlay
+
+
+@st.composite
+def _birefringence(draw):
+    shape = (2,) + tuple(
+        draw(st.lists(st.integers(1, 16), min_size=2, max_size=4))
+    )
+    dtype = draw(npst.floating_dtypes(sizes=(32, 64)))
+    bit_width = dtype.itemsize * 8
+    retardance = draw(
+        npst.arrays(
+            dtype,
+            shape=shape,
+            elements=st.floats(
+                min_value=0,
+                max_value=50,
+                exclude_min=True,
+                width=bit_width,
+            ),
+        )
+    )
+    orientation = draw(
+        npst.arrays(
+            dtype,
+            shape=shape,
+            elements=st.floats(
+                min_value=0,
+                max_value=dtype.type(np.pi),
+                exclude_min=True,
+                exclude_max=True,
+                width=bit_width,
+            ),
+        )
+    )
+    return retardance, orientation
+
+
+@given(briefringence=_birefringence(), jch=st.booleans())
+def test_ret_ori_overlay(briefringence: tuple[NDArray, NDArray], jch: bool):
+    """Test recOrder.io.utils.ret_ori_overlay()"""
+    retardance, orientation = briefringence
+    retardance_copy = retardance.copy()
+    orientation_copy = orientation.copy()
+    cmap = "JCh" if jch else "HSV"
+    overlay = ret_ori_overlay(
+        retardance,
+        orientation,
+        ret_max=np.percentile(retardance, 99),
+        cmap=cmap,
+    )
+    # check that the function did not mutate input data
+    assert_equal(retardance, retardance_copy)
+    assert_equal(orientation, orientation_copy)
+    # check output properties
+    # output contains NaN, pending further investigation
+    # assert overlay.min() >= 0
+    # assert overlay.max() <= 1
+    assert overlay.shape == retardance.shape + (3,)