diff --git a/linumpy/microscope/oct.py b/linumpy/microscope/oct.py
index 72e01166..2cda45f3 100644
--- a/linumpy/microscope/oct.py
+++ b/linumpy/microscope/oct.py
@@ -102,6 +102,36 @@ def load_image(self, crop: bool = True, galvo_shift: Union[bool, int] = False, c
             vol = vol[:, 0:n_alines, 0:n_bscans]
         return vol
 
+    def load_fringe(self, crop: bool = True):
+        # Create numpy array
+        # n_avg = self.info['n_repeat']  # TODO: use the number of averages when loading the data
+        n_alines = self.info['nx']
+        n_bscans = self.info['ny']
+        n_extra = self.info['n_extra']
+        n_alines_per_bscan = n_alines + n_extra
+        n_z = 2048
+
+        # Load the fringe
+        files = list(self.directory.rglob("fringe_*.bin"))
+        files.sort()
+        vol = None
+        for file in files:
+            with open(file, "rb") as f:
+                foo = np.fromfile(f, dtype=np.uint16)
+            n_frames = int(len(foo) / (n_alines_per_bscan * n_z))
+            foo = np.reshape(foo, (n_z, n_alines_per_bscan, n_frames), order='F')
+            if vol is None:
+                vol = foo
+            else:
+                vol = np.concatenate((vol, foo), axis=2)
+
+        # Crop the volume
+        if crop:
+            vol = vol[:, 0:n_alines, 0:n_bscans]
+
+
+        return vol
+
     def detect_galvo_shift(self, vol: np.ndarray = None) -> int:
         """Detect the galvo shift necessary to place the galvo return at the end of the bscans stack"""
         if vol is None:
diff --git a/scripts/dev/dev_convert_data_for_bouma.py b/scripts/dev/dev_convert_data_for_bouma.py
new file mode 100644
index 00000000..3959052f
--- /dev/null
+++ b/scripts/dev/dev_convert_data_for_bouma.py
@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""Prepare the data for the OCT k-space linearization and dispersion calibration."""
+
+import argparse
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.io import savemat
+from skimage.filters import threshold_li
+from tqdm.auto import tqdm
+
+from linumpy.microscope.oct import OCT
+
+# Parameters
+# directory = Path(r"D:\joel\2024-11-15-SOCT-UQAM-Calibration-Dispersion-Air-Mirror-woWindow-10x-higherContrast")
+
+def _build_arg_parser():
+    p = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawTextHelpFormatter)
+    p.add_argument("directory",
+                   help="Full path to a directory containing multiple OCT acquisition of a reference mirror at multiple depths.")
+    p.add_argument("basename", default="mirror_z",
+                   help="OCT base name, used to detect the files to process (default=%(default)s)")
+
+    return p
+def main():
+    # Parse arguments
+    p = _build_arg_parser()
+    args = p.parse_args()
+
+    # Parameters
+    directory = Path(args.directory)
+    basename = args.basename
+    spectrum_margin = 100
+
+    fringes_files = list(directory.glob(f"{basename}*"))
+    print(f"There are {len(fringes_files)} fringes files to process")
+
+    # Detect the file number
+    fringes_files.sort()
+    fringes_ids = [int(f.stem.split("_")[-1]) for f in fringes_files]
+
+    # Load the first fringe and determine the size of the data
+    oct = OCT(fringes_files[0])
+    galvo_return = oct.info['n_extra']
+    fringe = oct.load_fringe()
+    n_samples = fringe.shape[0]
+    n_alines = fringe.shape[1]
+    n_depths = len(fringes_files)
+    data = np.zeros((n_samples, n_alines, n_depths), dtype=fringe.dtype)
+    print(f"The dataset will have a shape of {data.shape}")
+
+    # Load all the data
+    for i, file in tqdm(zip(fringes_ids, fringes_files), total=n_depths):
+        oct = OCT(file)
+        fringe = oct.load_fringe()
+        data[..., i] = fringe[..., 0]
+
+    # Removing galvo return
+    data = data[:, 0:-galvo_return, ...]
+
+    # Detect the min/max wavelentgh with signal, and propose cropping
+    fringe_avg = data.mean(axis=(1, 2))
+    mask = fringe_avg > threshold_li(fringe_avg)
+    low_pt_otsu = max(np.where(mask)[0][0] - spectrum_margin, 0)
+    high_pt_otsu = min(np.where(mask)[0][-1] + spectrum_margin, n_samples)
+
+    figure_filename = directory / "spectrum.png"
+    plt.title(f"Detected spectrum from {low_pt_otsu} to {high_pt_otsu} ({high_pt_otsu - low_pt_otsu} samples)")
+    plt.plot(data.mean(axis=(1, 2)), label="Average fringe")
+    plt.axvspan(0, low_pt_otsu, color='r', alpha=0.5, label="Cropping region")
+    plt.axvspan(high_pt_otsu, n_samples, color='r', alpha=0.5)
+    plt.xlabel("Wavelength (px)")
+    plt.ylabel("Intensity")
+    plt.legend()
+    plt.savefig(figure_filename, dpi=300, pad_inches=0, bbox_inches='tight')
+    plt.show()
+
+    # Removing the background for each bscan
+    data_wo_bg = np.zeros_like(data, dtype=np.float64)
+    for i in range(data.shape[2]):
+        data_wo_bg[..., i] = data[..., i] - data[..., i].mean(axis=1, keepdims=True)
+
+    # Perform a quick reconstruction
+    window = np.hanning(data_wo_bg.shape[0]).reshape((data.shape[0], 1, 1))
+    bscans = np.abs(np.log(np.fft.fft(data_wo_bg * window, axis=0)))
+    bscans = bscans[0: n_samples // 2, :]
+
+    # Save the figure.
+    figure_filename = directory / "mirror_woCalibration.png"
+    plt.imshow(bscans.mean(axis=1), aspect="auto")
+    plt.xlabel("Mirror depth")
+    plt.ylabel("B-Scan Depth")
+    plt.title("Tomogram")
+    plt.savefig(figure_filename, dpi=300, bbox_inches='tight', pad_inches=0)
+    plt.show()
+
+    # Export the fringes as a matlab file
+    output = {'fringes': {'ch1': data_wo_bg.astype(np.float32)}, 'omitBscansVec': []}
+    output_file = directory / 'fringes.mat'
+    savemat(output_file, output)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/dev/dev_export_calibration_bouma.py b/scripts/dev/dev_export_calibration_bouma.py
new file mode 100644
index 00000000..31dc3ee9
--- /dev/null
+++ b/scripts/dev/dev_export_calibration_bouma.py
@@ -0,0 +1,85 @@
+# Exports the calibration files required by our OCTto be used with PolyMtl/UQAM OCTs
+from datetime import datetime
+from pathlib import Path
+
+from matplotlib import pyplot as plt
+from scipy.io import loadmat
+import numpy as np
+import json
+
+directory = Path("C://Users//Joël//OneDrive - UQAM//Documents//MATLAB//RobustMapping//Output//Results")
+output_directory = Path("C://Users//Public//Documents")
+result_file = directory / "UQAM-SOCT-20241115-Simplex-cropped.mat"
+
+name = "2024-11-29-air-10x"
+kmin = 852
+kmax = 1620
+n_samples = kmax - kmin
+oversampling_factor = 8
+n_samples_kspace = oversampling_factor * n_samples
+n_samples_output = 1024
+
+# Load the calibration optimization results
+calibration = loadmat(result_file)
+
+# Process the k-space mapping
+kspace_mapping = calibration["map"].squeeze()
+kspace_mapping = kspace_mapping / max(kspace_mapping) * n_samples_kspace # Adapting to the oversampling factor
+kspace_interp_matrix = np.zeros((n_samples//2, n_samples_kspace))
+for i in range(n_samples_kspace):
+    x = np.arange(n_samples_kspace)
+    y = np.zeros((n_samples_kspace,))
+    y[i] = 1
+    x_new = kspace_mapping
+
+    kspace_interp_matrix[:, i] = np.interp(x_new, x, y)
+
+# Process the dispersion
+dispersion = calibration["dispFit"].squeeze()
+x = np.linspace(-1, 1, len(dispersion))
+x_new = np.linspace(-1, 1, n_samples//2)
+dispersion = np.interp(x_new, x, dispersion)
+
+# Prepare the window
+window = np.hanning(n_samples//2)
+
+# Prepare the output
+figure_filename = output_directory / f"calibration_graphs_{name}.png"
+fig, axes = plt.subplots(nrows=1, ncols=3, figsize=[12,4])
+axes[0].plot(kspace_mapping)
+axes[0].set_title("K-Space Mapping")
+axes[1].plot(dispersion)
+axes[1].set_title("Dispersion")
+axes[2].plot(window)
+axes[2].set_title("Window")
+plt.savefig(figure_filename, dpi=300, pad_inches=0, bbox_inches='tight')
+
+
+# Write the calibration files
+window_filename = output_directory / f"filter_{name}.dat"
+with open(window_filename, "wb") as f:
+    window.astype(np.float64).tofile(f)
+
+dispersion_filename = output_directory / f"phase_{name}.dat"
+with open(dispersion_filename, "wb") as f:
+    dispersion.astype(np.float64).tofile(f)
+
+kspace_interp_matrix_filename = output_directory / f"kspace_interp_matrix_{name}.dat"
+with open(kspace_interp_matrix_filename, "wb") as f:
+    kspace_interp_matrix.astype(np.float64).tofile(f)
+
+config = {}
+config["name"] = name
+config["kmin"] = kmin
+config["kmax"] = kmax
+config["oversampling_factor"] = oversampling_factor
+config["n_samples_kspace"] = n_samples_kspace
+config["n_samples_output"] = n_samples_output
+config["window_filename"] = str(window_filename)
+config["dispersion_filename"] = str(dispersion_filename)
+config["kspace_interp_matrix_filename"] = str(kspace_interp_matrix_filename)
+config["calib_graph_filename"] = str(figure_filename)
+config["date"] = datetime.now().strftime("%Y-%m-%d@%Hh%M:%Ss")
+config_filename = output_directory / f"config_{name}.json"
+with open(config_filename, "w") as f:
+    json.dump(config, f)
\ No newline at end of file
diff --git a/scripts/dev_reconstruct_oct_bouma.py b/scripts/dev_reconstruct_oct_bouma.py
new file mode 100644
index 00000000..da5d1862
--- /dev/null
+++ b/scripts/dev_reconstruct_oct_bouma.py
@@ -0,0 +1,184 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from linumpy.microscope.oct import OCT
+from pathlib import Path
+from scipy.interpolate import interp1d
+import nibabel as nib
+from tqdm.auto import tqdm
+from scipy.io import loadmat
+
+#directory = Path(r"D:\joel\2024-11-15-SOCT-UQAM-Calibration-Dispersion-Air-Mirror-woWindow-10x-higherContrast\mirror_z_015")
+directory = Path(r"D:\joel\2024-11-15-SOCT-UQAM-Calibration-Dispersion-Air-Mirror-woWindow-10x-higherContrast")
+#directory = Path(r"D://joel//2024-11-01-SOCT-UQAM-Finger")
+#fringe_file = directory / "fringe"
+fringe_file = directory / 'fringes.mat'
+output_file = directory / "tomo_widerSpectrum.nii"
+
+# Parameters
+suffix = "2024-11-29-air-10x"
+calibration_dir = Path(r"C://Users//Public//Documents")
+kspace_interpolation_file = calibration_dir / f"kspace_interp_matrix_{suffix}.dat"
+kspace_interp_file = calibration_dir / f"kspace_interp_matrix_{suffix}.dat"
+window_file = calibration_dir / f"filter_{suffix}.dat"
+dispersion_file = calibration_dir / f"phase_{suffix}.dat"
+
+replace_by_default = False
+
+fringe = loadmat(fringe_file)['fringes'][0, 0][0]
+# galvo_return = 40
+# spectrum_margin = 100
+#oct = OCT(fringe_file)
+
+# Load the data and calibration files
+#fringe = oct.load_fringe(crop=False).squeeze()
+n_samples = fringe.shape[0]
+n_alines = fringe.shape[1]
+n_bscans = fringe.shape[2]
+n_samples = 2048
+oversamplingFactor = 8
+k_min = 852
+k_max = 1620
+n_values = k_max - k_min
+oversampled_kspace = n_values * oversamplingFactor
+dc_margin = 0
+output_size = 1024
+print(f"The data shape is {fringe.shape}")
+
+# Load the calibration files
+window = np.fromfile(window_file, dtype=np.float64)
+dispersion = np.fromfile(dispersion_file, dtype=np.float64)
+kspace_sparse_interp = np.fromfile(kspace_interpolation_file, dtype=np.float64)
+kspace_sparse_interp = np.reshape(kspace_sparse_interp, (len(dispersion), len(kspace_sparse_interp) // len(dispersion)), order="C")
+
+map = np.fromfile(kspace_interp_file, dtype=np.float64)
+
+
+def reconstruct(fringe, show: bool = False):
+    # 1. Input fringes of shape N_Samples x N_Alines
+    if show:
+        plt.imshow(fringe)
+        plt.axhline(k_min, color='red', linestyle='--')
+        plt.axhline(k_max, color='red', linestyle='--')
+        plt.show()
+
+    # 2. Crop the wavelengths to remove empty spectral regions
+    fringe = fringe[k_min:k_max, :]
+    if show:
+        plt.imshow(fringe)
+        plt.title("2. Cropped fringe")
+        plt.show()
+
+    # 3. Compute (update) the background from the cropped fringe
+    background = fringe.mean(axis=1, keepdims=True)
+    if show:
+        plt.plot(background)
+        plt.title("3. Background")
+        plt.show()
+
+    # 4. Remove background
+    fringe = fringe - background
+    if show:
+        plt.imshow(fringe)
+        plt.title("4. Fringes w/o background")
+        plt.show()
+
+    # 5. Compte a preprocessing tomograph
+    preTom = np.fft.fft(fringe, axis=0)
+    if show:
+        plt.imshow(np.abs(np.log(preTom)))
+        plt.title("5. PreTom")
+        plt.show()
+
+    # 6. Mask the complex conjugate and (Optional) 7. Mask the DC (optional)
+    preTom[0:n_values // 2, :] = 0
+    if dc_margin > 0:
+        preTom[-dc_margin::, :] = 0
+    if show:
+        plt.imshow(np.abs(np.log(preTom + 1e-6)))
+        plt.title("6/7. Masked PreTom")
+        plt.show()
+
+    # 8. Zero-padding to oversample the k-space for the interpolation
+    padding_size = oversampled_kspace - preTom.shape[0]
+    preTom = np.pad(preTom, ((padding_size, 0), (0, 0)), mode='constant', constant_values=0)
+    if show:
+        plt.imshow(np.abs(np.log(preTom + 1e-6)))
+        plt.title("8. Padded PreTom")
+        plt.show()
+
+    # 9. Go back to k-space
+    fringe = np.fft.ifft(preTom, axis=0)
+    if show:
+        plt.imshow(np.abs(fringe))
+        plt.title("9. Oversampled Complex Fringes")
+        plt.show()
+
+    # 10. K-space linearization via sparse matrix interpolation
+    fringe_linear = kspace_sparse_interp @ fringe
+    # f = interp1d(np.linspace(0, n_samples, n_samples), fringe, axis=0)
+    # fringe_linear = f(map)
+    if show:
+        plt.imshow(np.abs(fringe_linear))
+        plt.title("10. Lineararized K-Space Fringe")
+        plt.show()
+
+    # 11. Compensate dispersion and apodization
+    complex_window = window * np.exp(-1j * dispersion)
+    x = np.linspace(0, 1, complex_window.shape[0])
+    y = complex_window.squeeze()
+    f = interp1d(x, y)
+    complex_window = f(np.linspace(0, 1, fringe_linear.shape[0]))
+    fringe_linear = fringe_linear * complex_window.reshape([*complex_window.shape, 1])
+
+    if show:
+        plt.imshow(np.abs(fringe_linear))
+        plt.title("11. Lineararized K-Space Fringe")
+        plt.show()
+
+    # 12. Zero-padding to the output size
+    #output_size = n_values // 2
+    padding_size = output_size - fringe_linear.shape[0]
+    fringe_linear = np.pad(fringe_linear, ((0, padding_size), (0, 0)))
+    if show:
+        plt.imshow(np.abs(fringe_linear))
+        plt.title("12. Lineararized K-Space Fringe, padded")
+        plt.show()
+
+    # 13. Circular shift
+    shift = kspace_sparse_interp.shape[0] // 2
+    fringe_linear = np.roll(fringe_linear, -shift, axis=0)
+    if show:
+        plt.imshow(np.abs(fringe_linear))
+        plt.title("13. Centered processed fringes")
+        plt.show()
+
+    # 14. Final FFT
+    tomogram = np.fft.fft(fringe_linear, axis=0)
+    if show:
+        plt.imshow(np.abs(np.log(tomogram + 1e-6)))
+        plt.title("14. Tomogram")
+        plt.show()
+
+    return tomogram
+
+
+n_samples = n_values // 2
+n_alines = fringe.shape[1]
+tomogram = np.zeros((output_size, n_alines, n_bscans), dtype=np.float32)
+for i in tqdm(range(n_bscans)):
+    tomo = reconstruct(fringe[:, :, i].squeeze(), show=False)
+    tomogram[:, :, i] = np.abs(tomo)
+tomogram = np.flipud(tomogram)
+
+figure_filename = directory / "mirror_calibrated.png"
+plt.imshow(np.log(tomogram.mean(axis=1)), aspect="auto")
+plt.xlabel("Mirror depth")
+plt.ylabel("B-Scan Depth")
+plt.title("Calibrated Tomograms")
+plt.savefig(figure_filename, dpi=300, bbox_inches='tight', pad_inches=0)
+plt.show()
+
+# Save the volume
+# affine = np.eye(4)
+# img = nib.Nifti1Image(tomogram, affine)
+# nib.save(img, output_file)