diff --git a/WrapImage/nifti_wrapper.py b/WrapImage/nifti_wrapper.py
index 2117e57d..ba2339c2 100644
--- a/WrapImage/nifti_wrapper.py
+++ b/WrapImage/nifti_wrapper.py
@@ -100,7 +100,7 @@ def loop_over_first_n_minus_1_dimensions(arr):
 
         # Pass additional arguments to the algorithm
 
-        fit = OsipiBase(algorithm=args.algorithm)
+        fit = OsipiBase(algorithm=args.algorithm, bvalues=bvals)
         f_image = []
         Dp_image = []
         D_image = []
@@ -109,7 +109,7 @@ def loop_over_first_n_minus_1_dimensions(arr):
         n = data.ndim
         total_iteration = np.prod(data.shape[:n-1])
         for idx, view in tqdm(loop_over_first_n_minus_1_dimensions(data), desc=f"{args.algorithm} is fitting", dynamic_ncols=True, total=total_iteration):
-            fit_result = fit.osipi_fit(view, bvals)
+            fit_result = fit.osipi_fit(view)
             f_image.append(fit_result["f"])
             Dp_image.append(fit_result["Dp"])
             D_image.append(fit_result["D"])
diff --git a/WrapImage/nifti_wrapper_kaapana.py b/WrapImage/nifti_wrapper_kaapana.py
index b4a8856f..f0fd96c8 100644
--- a/WrapImage/nifti_wrapper_kaapana.py
+++ b/WrapImage/nifti_wrapper_kaapana.py
@@ -127,7 +127,7 @@ def load_config():
         affine = np.array(affine_override).reshape(4, 4)
 
     # Initialize model
-    fit = OsipiBase(algorithm=algorithm)
+    fit = OsipiBase(algorithm=algorithm, bvalues=bvals)
 
     # Preallocate output arrays
     shape = data.shape[:data.ndim - 1]
@@ -142,7 +142,7 @@ def load_config():
         loop_over_first_n_minus_1_dimensions(data),
         desc="Fitting IVIM model", dynamic_ncols=True, total=total_iteration
     ):
-        fit_result = fit.osipi_fit(view, bvals)
+        fit_result = fit.osipi_fit(view)
         f_image[idx] = fit_result["f"]
         Dp_image[idx] = fit_result["Dp"]
         D_image[idx] = fit_result["D"]
diff --git a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
index 6ee3735c..9c49e859 100644
--- a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
+++ b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
@@ -66,12 +66,11 @@ def algorithm(self,dwi_arr, bval_arr, LB0, UB0, x0in):
         (f_arr, D_arr, Dx_arr, s0_arr, fitted_dwi_arr, RSS, rms_val, chi, AIC, BIC, R_sq) = results
         return D_arr/1000, f_arr, Dx_arr/1000, s0_arr
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -81,7 +80,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         LB = np.array(bounds[0])[[1,0,2,3]]
         UB = np.array(bounds[1])[[1,0,2,3]]
 
diff --git a/src/standardized/DT_IIITN_WLS.py b/src/standardized/DT_IIITN_WLS.py
index 97f1fb9a..675ef5d2 100644
--- a/src/standardized/DT_IIITN_WLS.py
+++ b/src/standardized/DT_IIITN_WLS.py
@@ -66,17 +66,16 @@ def __init__(self, bvalues=None, thresholds=None,
         )
         self.method = method.upper()
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit using the selected method (WLS or RLM).
 
         Args:
             signals (array-like): Signal intensities at each b-value.
-            bvalues (array-like, optional): b-values for the signals.
 
         Returns:
             dict: Dictionary with keys "D", "f", "Dp".
         """
-        bvalues = np.array(bvalues)
+        bvalues = np.array(self.bvalues)
 
         # Use threshold as cutoff if available
         cutoff = 200
diff --git a/src/standardized/ETP_SRI_LinearFitting.py b/src/standardized/ETP_SRI_LinearFitting.py
index f963c014..873396d7 100644
--- a/src/standardized/ETP_SRI_LinearFitting.py
+++ b/src/standardized/ETP_SRI_LinearFitting.py
@@ -59,20 +59,18 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         # Check the inputs
         
     
-    def ivim_fit(self, signals, bvalues=None, linear_fit_option=False, **kwargs):
+    def ivim_fit(self, signals, linear_fit_option=False, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
             linear_fit_option (bool, optional): This fit has an option to only run a linear fit. Defaults to False.
 
         Returns:
             _type_: _description_
         """
         signals[signals<0.0000001]=0.0000001
-        if bvalues is None:
-            bvalues = self.bvalues
+        bvalues = self.bvalues
         
         if self.thresholds is None:
             ETP_object = LinearFit()
diff --git a/src/standardized/IAR_LU_biexp.py b/src/standardized/IAR_LU_biexp.py
index e362bf1a..b7928f71 100644
--- a/src/standardized/IAR_LU_biexp.py
+++ b/src/standardized/IAR_LU_biexp.py
@@ -71,12 +71,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -88,10 +87,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
@@ -110,12 +106,11 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         return results
 
 
-    def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
+    def ivim_fit_full_volume(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -125,16 +120,14 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
                        [self.bounds["S0"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["D"][1]]]
         initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
             gtab = gradient_table(bvalues, bvecs=bvec, b0_threshold=0)
             
             self.IAR_algorithm = IvimModelBiExp(gtab, bounds=bounds, initial_guess=initial_guess)
+        bvalues = np.asarray(self.bvalues)
         b0_index = np.where(bvalues == 0)[0][0]
         mask = signals[...,b0_index]>0
         fit_results = self.IAR_algorithm.fit(signals, mask=mask)
diff --git a/src/standardized/IAR_LU_modified_mix.py b/src/standardized/IAR_LU_modified_mix.py
index 48aef405..473464d4 100644
--- a/src/standardized/IAR_LU_modified_mix.py
+++ b/src/standardized/IAR_LU_modified_mix.py
@@ -68,12 +68,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -83,10 +82,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
                   [self.bounds["f"][1], self.bounds["Dp"][1]*1000, self.bounds["D"][1]*1000]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
diff --git a/src/standardized/IAR_LU_modified_topopro.py b/src/standardized/IAR_LU_modified_topopro.py
index d86b6779..d4f6509d 100644
--- a/src/standardized/IAR_LU_modified_topopro.py
+++ b/src/standardized/IAR_LU_modified_topopro.py
@@ -72,12 +72,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -86,10 +85,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
                   [self.bounds["f"][1], self.bounds["Dp"][1]*1000, self.bounds["D"][1]*1000]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
diff --git a/src/standardized/IAR_LU_segmented_2step.py b/src/standardized/IAR_LU_segmented_2step.py
index dbbec2c9..9831efa0 100644
--- a/src/standardized/IAR_LU_segmented_2step.py
+++ b/src/standardized/IAR_LU_segmented_2step.py
@@ -72,12 +72,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, thresholds=None, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -90,10 +89,7 @@ def ivim_fit(self, signals, bvalues, thresholds=None, **kwargs):
         initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
diff --git a/src/standardized/IAR_LU_segmented_3step.py b/src/standardized/IAR_LU_segmented_3step.py
index 5a03bee6..df16f6b5 100644
--- a/src/standardized/IAR_LU_segmented_3step.py
+++ b/src/standardized/IAR_LU_segmented_3step.py
@@ -75,12 +75,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -93,10 +92,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
diff --git a/src/standardized/IAR_LU_subtracted.py b/src/standardized/IAR_LU_subtracted.py
index c2b6fd3b..890b8be1 100644
--- a/src/standardized/IAR_LU_subtracted.py
+++ b/src/standardized/IAR_LU_subtracted.py
@@ -72,12 +72,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             self.IAR_algorithm = None
         
     
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -90,10 +89,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         
         if self.IAR_algorithm is None:
-            if bvalues is None:
-                bvalues = self.bvalues
-            else:
-                bvalues = np.asarray(bvalues)
+            bvalues = np.asarray(self.bvalues)
             
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
diff --git a/src/standardized/IVIM_NEToptim.py b/src/standardized/IVIM_NEToptim.py
index e7c77193..0595e3bb 100644
--- a/src/standardized/IVIM_NEToptim.py
+++ b/src/standardized/IVIM_NEToptim.py
@@ -89,18 +89,15 @@ def initialize(self, bounds, initial_guess, fitS0, traindata, SNR, n):
         self.algorithm =lambda data: deep.predict_IVIM(data, self.bvalues, self.net, self.arg)
 
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
-        if not np.array_equal(bvalues, self.bvalues):
-            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
 
         paramsNN = deep.predict_IVIM(signals, self.bvalues, self.net, self.arg)
 
@@ -112,18 +109,15 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         return results
 
 
-    def ivim_fit_full_volume(self, signals, bvalues, retrain_on_input_data=False, **kwargs):
+    def ivim_fit_full_volume(self, signals, retrain_on_input_data=False, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
-        if not np.array_equal(bvalues, self.bvalues):
-            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
         minimum_bvalue = np.min(self.bvalues) # We normalize the signal to the minimum bvalue. Should be 0 or very close to 0.
         b0_indices = np.where(self.bvalues == minimum_bvalue)[0]
         normalization_factor = np.mean(signals[..., b0_indices],axis=-1)
diff --git a/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py b/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
index 2de8cbf7..5a161556 100644
--- a/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
+++ b/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
@@ -73,12 +73,11 @@ def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None,
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2])
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
+    def ivim_fit(self, signals, initial_guess=None, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -88,7 +87,7 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
 
         initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
 
         epsilon = 0.000001
         fit_results = fit_segmented(bvalues, signals, bounds=bounds, cutoff=self.thresholds, p0=initial_guess)
@@ -105,11 +104,10 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
 
         return results
 
-    def ivim_fit_full_volume(self, signals, bvalues, njobs=4, **kwargs):
+    def ivim_fit_full_volume(self, signals, njobs=4, **kwargs):
         """Perform the IVIM fit
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
         Returns:
             _type_: _description_
         """
@@ -122,6 +120,7 @@ def ivim_fit_full_volume(self, signals, bvalues, njobs=4, **kwargs):
         # Get index of b=0
         shape=np.shape(signals)
 
+        bvalues=np.array(self.bvalues)
         b0_index = np.where(bvalues == 0)[0][0]
         # Mask of voxels where signal at b=0 >= 0.5
         valid_mask = signals[..., b0_index] >= 0
@@ -133,7 +132,7 @@ def ivim_fit_full_volume(self, signals, bvalues, njobs=4, **kwargs):
         normalization_factor = np.mean(signals[..., b0_indices],axis=-1)
         signals = signals / np.repeat(normalization_factor[...,np.newaxis],np.shape(signals)[-1],-1)
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
 
         epsilon = 0.000001
         fit_results = np.array(fit_segmented_array(bvalues, signals, bounds=bounds, cutoff=self.thresholds, p0=initial_guess))
diff --git a/src/standardized/OGC_AmsterdamUMC_biexp.py b/src/standardized/OGC_AmsterdamUMC_biexp.py
index 80e5a8e0..0ec3f797 100644
--- a/src/standardized/OGC_AmsterdamUMC_biexp.py
+++ b/src/standardized/OGC_AmsterdamUMC_biexp.py
@@ -51,12 +51,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         self.use_initial_guess = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
         self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -66,7 +65,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         fit_results = self.OGC_algorithm(bvalues, signals, p0=initial_guess, bounds=bounds, fitS0=self.fitS0)
 
         results = {}
diff --git a/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py b/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
index 6d687feb..4a8f20d1 100644
--- a/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
+++ b/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
@@ -58,12 +58,11 @@ def initialize(self, thresholds):
         else:
             self.thresholds = thresholds
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -73,7 +72,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
         initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
-        bvalues=np.array(bvalues)
+        bvalues = np.array(self.bvalues)
         fit_results = self.OGC_algorithm(bvalues, signals, bounds=bounds, cutoff=self.thresholds, p0=initial_guess)
 
         results = {}
diff --git a/src/standardized/OJ_GU_bayesMATLAB.py b/src/standardized/OJ_GU_bayesMATLAB.py
index 72ccaa8f..c20ef75a 100644
--- a/src/standardized/OJ_GU_bayesMATLAB.py
+++ b/src/standardized/OJ_GU_bayesMATLAB.py
@@ -68,12 +68,11 @@ def algorithm(self, Y, b, lim, blim, initial_guess):
         out = self.eng.IVIM_bayes(Y, f, D, Dstar, S0, b, lim, nargout=1)
         return out['D']['mode'], out['f']['mode'], out['Dstar']['mode'], out['S0']['mode']
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -84,7 +83,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
         fit_results = self.algorithm(np.array(signals)[:,np.newaxis], 
-                                     np.array(bvalues), 
+                                     np.array(self.bvalues), 
                                      np.array(bounds)[:,[1,0,2,3]], 
                                      self.thresholds,
                                      initial_guess)
diff --git a/src/standardized/OJ_GU_seg.py b/src/standardized/OJ_GU_seg.py
index 9b8cbf2f..24b584d9 100644
--- a/src/standardized/OJ_GU_seg.py
+++ b/src/standardized/OJ_GU_seg.py
@@ -49,21 +49,17 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         # Initialize the algorithm
         
     
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
         
-        if bvalues is None:
-            bvalues = self.bvalues
-        else:
-            bvalues = np.asarray(bvalues)
+        bvalues = np.asarray(self.bvalues)
 
         if self.thresholds is None:
             bthr = 200
diff --git a/src/standardized/OJ_GU_segMATLAB.py b/src/standardized/OJ_GU_segMATLAB.py
index e2d68125..7ee93927 100644
--- a/src/standardized/OJ_GU_segMATLAB.py
+++ b/src/standardized/OJ_GU_segMATLAB.py
@@ -65,12 +65,11 @@ def algorithm(self, Y, b, lim, blim):
         return pars['D'], pars['f'], pars['Dstar'], pars['S0']
 
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -79,7 +78,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
                   [self.bounds["D"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["S0"][1]])
 
         fit_results = self.algorithm(np.array(signals)[:,np.newaxis], 
-                                     np.array(bvalues), 
+                                     np.array(self.bvalues), 
                                      np.array(bounds)[:,[0,3,1,2]], 
                                      self.thresholds)
 
diff --git a/src/standardized/PV_MUMC_biexp.py b/src/standardized/PV_MUMC_biexp.py
index fa854cda..6cac6c28 100644
--- a/src/standardized/PV_MUMC_biexp.py
+++ b/src/standardized/PV_MUMC_biexp.py
@@ -44,26 +44,16 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         self.use_initial_guess = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
         
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             dict: Fitted IVIM parameters f, Dp (D*), and D.
         """
-        # --- bvalues resolution ---
-        # Edge case: bvalues not passed here → fall back to the ones set at __init__ time
-        if bvalues is None:
-            if self.bvalues is None:
-                raise ValueError(
-                    "PV_MUMC_biexp: bvalues must be provided either at initialization or at fit time."
-                )
-            bvalues = self.bvalues
-        else:
-            bvalues = np.asarray(bvalues)
+        bvalues = np.asarray(self.bvalues)
 
         # --- Bounds resolution ---
         # self.bounds is always a dict (OsipiBase force_default_settings=True).
diff --git a/src/standardized/PvH_KB_NKI_IVIMfit.py b/src/standardized/PvH_KB_NKI_IVIMfit.py
index a574a94e..22bee5cc 100644
--- a/src/standardized/PvH_KB_NKI_IVIMfit.py
+++ b/src/standardized/PvH_KB_NKI_IVIMfit.py
@@ -52,18 +52,17 @@ def __init__(self, bvalues=None, thresholds=None,bounds=None,initial_guess=None)
         self.use_initial_guess = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
 
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
         #bvalues = np.array(bvalues)
-        bvalues = bvalues.tolist() #NKI code expects a list instead of nparray
+        bvalues = np.asarray(self.bvalues).tolist() #NKI code expects a list instead of nparray
         # reshape signal as the NKI code expects a 4D array
         signals[signals<0.00001]=0.00001
         signals = np.reshape(signals, (1, 1, 1, len(signals)))  # assuming that in this test the signals are always single voxel
diff --git a/src/standardized/Super_IVIM_DC.py b/src/standardized/Super_IVIM_DC.py
index cc853e13..502abcea 100644
--- a/src/standardized/Super_IVIM_DC.py
+++ b/src/standardized/Super_IVIM_DC.py
@@ -88,18 +88,15 @@ def initialize(self, bounds, initial_guess, fitS0, SNR, working_dir=os.getcwd(),
         )
 
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             results: a dictionary containing "d", "f", and "Dp".
         """
-        if not np.array_equal(bvalues, self.bvalues):
-            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
 
         Dp, Dt, f, S0_superivimdc = infer_from_signal(
             signal=signals,
@@ -115,18 +112,15 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         return results
 
 
-    def ivim_fit_full_volume(self, signals, bvalues, retrain_on_input_data=False, **kwargs):
+    def ivim_fit_full_volume(self, signals, retrain_on_input_data=False, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
-        if not np.array_equal(bvalues, self.bvalues):
-            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
 
         nanmask = np.any(np.isnan(signals),axis=-1)
         signals,shape = self.reshape_to_voxelwise(signals)
diff --git a/src/standardized/TCML_TechnionIIT_SLS.py b/src/standardized/TCML_TechnionIIT_SLS.py
index 991b2a5d..e433b685 100644
--- a/src/standardized/TCML_TechnionIIT_SLS.py
+++ b/src/standardized/TCML_TechnionIIT_SLS.py
@@ -59,18 +59,17 @@ def initialize(self, bounds, fitS0,thresholds):
         else:
             self.thresholds = thresholds
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         bounds = ([self.bounds["D"][0], self.bounds["Dp"][0], self.bounds["f"][0], self.bounds["S0"][0]],
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         signals[signals<0]=0
diff --git a/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py b/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py
index 212a2cd0..5004e27a 100644
--- a/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py
+++ b/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py
@@ -59,18 +59,17 @@ def initialize(self, bounds, initial_guess, fitS0):
             self.use_initial_guess = {"f": True, "Dp": True, "D": True}
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         bounds = ([self.bounds["D"][0], self.bounds["Dp"][0], self.bounds["f"][0], self.bounds["S0"][0]],
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         initial_guess = [self.initial_guess["D"], self.initial_guess["Dp"], self.initial_guess["f"], self.initial_guess["S0"]]
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py b/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
index 7c9c81ef..c579e3d3 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
@@ -58,18 +58,17 @@ def initialize(self, bounds, fitS0, thresholds):
         self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
         self.use_initial_guess = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
         signals[signals<0]=0
-        bvalues=np.array(bvalues)
+        bvalues = np.array(self.bvalues)
         bounds = ([self.bounds["D"][0], self.bounds["Dp"][0], self.bounds["f"][0], self.bounds["S0"][0]],
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py b/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
index 065fb92b..a3528a08 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
@@ -60,18 +60,17 @@ def initialize(self, bounds, fitS0,thresholds):
         self.fitS0=fitS0
         self.use_initial_guess = {"f": False, "Dp": False, "D": False}
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
         signals[signals<0]=0
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         bounds = ([self.bounds["D"][0], self.bounds["Dp"][0], self.bounds["f"][0], self.bounds["S0"][0]],
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         fit_results = self.fit_least_squares(np.array(signals)[:,np.newaxis],bvalues, bounds, min_bval_high=self.thresholds)
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py b/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
index eeecd661..09a06fa5 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
@@ -62,18 +62,17 @@ def initialize(self, bounds, fitS0, thresholds):
             self.thresholds = thresholds
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
         signals[signals<0]=0
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         bounds = ([self.bounds["D"][0], self.bounds["Dp"][0], self.bounds["f"][0], self.bounds["S0"][0]],
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         fit_results = self.fit_least_squares(np.array(signals)[:,np.newaxis],bvalues, bounds,min_bval_high=self.thresholds)
diff --git a/src/standardized/TCML_TechnionIIT_lsqlm.py b/src/standardized/TCML_TechnionIIT_lsqlm.py
index 429fabcf..a8a57f1a 100644
--- a/src/standardized/TCML_TechnionIIT_lsqlm.py
+++ b/src/standardized/TCML_TechnionIIT_lsqlm.py
@@ -56,18 +56,17 @@ def initialize(self, bounds, initial_guess, fitS0):
             self.use_initial_guess = {"f": True, "Dp": True, "D": True}
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
         """
 
-        bvalues=np.array(bvalues)
+        bvalues=np.array(self.bvalues)
         initial_guess = [self.initial_guess["D"], self.initial_guess["Dp"], self.initial_guess["f"], self.initial_guess["S0"]]
         fit_results = self.fit_least_squares(bvalues, np.array(signals)[:,np.newaxis], initial_guess)
 
diff --git a/src/standardized/TCML_TechnionIIT_lsqtrf.py b/src/standardized/TCML_TechnionIIT_lsqtrf.py
index 64c0b6b3..fb0ccf8d 100644
--- a/src/standardized/TCML_TechnionIIT_lsqtrf.py
+++ b/src/standardized/TCML_TechnionIIT_lsqtrf.py
@@ -60,12 +60,11 @@ def initialize(self, bounds, initial_guess, fitS0):
 
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues, **kwargs):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -74,7 +73,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
                        [self.bounds["D"][1], self.bounds["Dp"][1], self.bounds["f"][1], self.bounds["S0"][1]])
         initial_guess = [self.initial_guess["D"], self.initial_guess["Dp"], self.initial_guess["f"], self.initial_guess["S0"]]
 
-        bvalues=np.array(bvalues)
+        bvalues = np.array(self.bvalues)
         fit_results = self.fit_least_squares(bvalues, np.array(signals)[:,np.newaxis], bounds,initial_guess)
 
         def get_scalar(val):
diff --git a/src/standardized/TF_reference_IVIMfit.py b/src/standardized/TF_reference_IVIMfit.py
index f72062d5..6114578a 100644
--- a/src/standardized/TF_reference_IVIMfit.py
+++ b/src/standardized/TF_reference_IVIMfit.py
@@ -57,12 +57,11 @@ def initialize(self, bounds, thresholds):
             self.thresholds = thresholds
 
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, **kwargs):
         """Perform the IVIM fit
 
         Args:
             signals (array-like)
-            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
 
         Returns:
             _type_: _description_
@@ -70,7 +69,7 @@ def ivim_fit(self, signals, bvalues=None):
         bounds = ([self.bounds["D"][0], self.bounds["f"][0], self.bounds["Dp"][0], self.bounds["S0"][0]],
                   [self.bounds["D"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["S0"][1]])
 
-        bvalues = np.array(bvalues)
+        bvalues = np.array(self.bvalues)
 
         if np.any(signals < 0):
             signals = np.clip(signals,0.01, None)
diff --git a/src/wrappers/OsipiBase.py b/src/wrappers/OsipiBase.py
index 3555dde7..2c775c02 100644
--- a/src/wrappers/OsipiBase.py
+++ b/src/wrappers/OsipiBase.py
@@ -80,10 +80,10 @@ class OsipiBase:
     osipi_initiate_algorithm(algorithm, **kwargs)
         Dynamically replace the current instance with the specified
         algorithm subclass.
-    osipi_fit(data, bvalues=None, njobs=1, **kwargs)
+    osipi_fit(data, njobs=1, **kwargs)
         Voxel-wise IVIM fitting with optional parallel processing and
         automatic signal normalization.
-    osipi_fit_full_volume(data, bvalues=None, **kwargs)
+    osipi_fit_full_volume(data, **kwargs)
         Full-volume fitting for algorithms that support it.
     osipi_print_requirements()
         Display algorithm requirements such as needed b-values or bounds.
@@ -91,7 +91,7 @@ class OsipiBase:
         Query acceptable input dimensionalities.
     osipi_check_required_*()
         Validate that provided inputs meet algorithm requirements.
-    osipi_simple_bias_and_RMSE_test(SNR, bvalues, f, Dstar, D, noise_realizations)
+    osipi_simple_bias_and_RMSE_test(SNR, f, Dstar, D, noise_realizations)
         Monte-Carlo bias/RMSE evaluation of the current fitting method.
     D_and_Ds_swap(results)
         Ensure consistency of D and D* estimates by swapping if necessary.
@@ -239,8 +239,7 @@ def initialize(**kwargs):
         """Placeholder for subclass initialization"""
         pass
 
-    #def osipi_fit(self, data=None, bvalues=None, thresholds=None, bounds=None, initial_guess=None, **kwargs):
-    def osipi_fit(self, data, bvalues=None, njobs=1, **kwargs):
+    def osipi_fit(self, data, njobs=1, **kwargs):
         """
         Fit multi-b-value diffusion MRI data using the IVIM model.
 
@@ -253,9 +252,6 @@ def osipi_fit(self, data, bvalues=None, njobs=1, **kwargs):
         data : np.ndarray
             Multi-dimensional array containing the signal intensities. The last dimension must correspond
             to the b-values (diffusion weightings).
-        bvalues : array-like, optional
-            Array of b-values corresponding to the last dimension of `data`. If not provided, the method
-            uses `self.bvalues` set during object initialization.
         njobs : int, optional, default=1
             Number of parallel jobs to use for voxel-wise fitting. If `njobs` > 1, the fitting will be
             distributed across multiple processes. -1 will use all available cpus
@@ -279,40 +275,18 @@ def osipi_fit(self, data, bvalues=None, njobs=1, **kwargs):
 
         Example
         -------
-        >>> results = instance.osipi_fit(data, bvalues=[0, 50, 200, 800], njobs=4)
+        >>> results = instance.osipi_fit(data, njobs=4)
         >>> f_map = results['f']
         >>> D_map = results['D']
         """
         
-        # We should first check whether the attributes in the __init__ are not None
-        # Then check if they are input here, if they are, these should overwrite the attributes
-        use_bvalues = bvalues if bvalues is not None else self.bvalues
-        #use_thresholds = thresholds if self.thresholds is None else self.thresholds
-        #use_bounds = bounds if self.bounds is None else self.bounds
-        #use_initial_guess = initial_guess if self.initial_guess is None else self.initial_guess
-
-        # Make sure we don't make arrays of None's
-        if use_bvalues is not None: use_bvalues = np.asarray(use_bvalues) 
-        #if use_thresholds is not None: use_thresholds = np.asarray(use_thresholds) 
-        #if use_bounds is not None: use_bounds = np.asarray(use_bounds) 
-        #if use_initial_guess is not None: use_initial_guess = np.asarray(use_initial_guess) 
-        #kwargs["bvalues"] = use_bvalues
-        
-        #args = [data, use_bvalues, use_thresholds]
-        #if self.required_bounds or self.required_bounds_optional:
-            #args.append(use_bounds)
-        #if self.required_initial_guess or self.required_initial_guess_optional:
-            #args.append(use_initial_guess)
-            
-        # Run a check_requirements method that makes sure that these inputs fulfil the requirements
-        
-        
-        # Then we check which inputs are required for the algorithm using the requirements attributes in the template 
-        
-        # Then we pass everything into the ivim_fit method which performs the fit according to the algorithm
-        
-        #args = [data, use_bvalues, use_initial_guess, use_bounds, use_thresholds]
-        #args = [arg for arg in args if arg is not None]
+        # b-values must be provided at initialization, not at fit-time.
+        if self.bvalues is None:
+            raise ValueError(
+                "b-values must be provided at initialization (__init__), not at fit-time. "
+                "Pass bvalues to the constructor: OsipiBase(bvalues=..., algorithm=...)"
+            )
+        use_bvalues = np.asarray(self.bvalues)
 
         # Check if there is an attribute that defines the result dictionary keys
         if hasattr(self, "result_keys"):
@@ -349,7 +323,7 @@ def osipi_fit(self, data, bvalues=None, njobs=1, **kwargs):
             def parfun(ijk):
                 single_voxel_data = np.array(data[ijk], copy=True)
                 if not np.isnan(single_voxel_data[0]):
-                    fit = self.ivim_fit(single_voxel_data, use_bvalues, **kwargs)
+                    fit = self.ivim_fit(single_voxel_data, **kwargs)
                 else:
                     fit={'D':0,'f':0,'Dp':0}
                 return ijk, fit
@@ -372,8 +346,7 @@ def parfun(ijk):
                 # Normalize array
                 single_voxel_data = data[ijk]
                 if not np.isnan(single_voxel_data[0]):
-                    args = [single_voxel_data, use_bvalues]
-                    fit = self.ivim_fit(*args, **kwargs)
+                    fit = self.ivim_fit(single_voxel_data, **kwargs)
                 else:
                     fit={'D':0,'f':0,'Dp':0}
                 for key in list(fit.keys()):
@@ -382,7 +355,7 @@ def parfun(ijk):
         return results
 
 
-    def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
+    def osipi_fit_full_volume(self, data, **kwargs):
         """
         Fit an entire volume of multi-b-value diffusion MRI data in a single call using the IVIM model.
 
@@ -394,9 +367,6 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
         data : np.ndarray
             2D (data x b-values), 3D (single slice), or 4D (multi-slice) diffusion-weighted imaging (DWI) data.
             The last dimension must correspond to the b-values.
-        bvalues : array-like, optional
-            Array of b-values corresponding to the last dimension of `data`. If not provided, the method
-            uses `self.bvalues` set during object initialization.
         **kwargs : dict, optional
             Additional keyword arguments to be passed to `ivim_fit_full_volume`.
 
@@ -414,7 +384,7 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
         Example
         -------
         # Standard usage:
-        results = instance.osipi_fit_full_volume(data, bvalues=[0, 50, 200, 800])
+        results = instance.osipi_fit_full_volume(data)
         f_map = results['f']
         D_map = results['D']
         Dp_map = results['Dp']
@@ -426,8 +396,12 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
         """
         
         try:
-            use_bvalues = bvalues if bvalues is not None else self.bvalues
-            if use_bvalues is not None: use_bvalues = np.asarray(use_bvalues) 
+            # b-values must be provided at initialization, not at fit-time.
+            if self.bvalues is None:
+                raise ValueError(
+                    "b-values must be provided at initialization (__init__), not at fit-time. "
+                    "Pass bvalues to the constructor: OsipiBase(bvalues=..., algorithm=...)"
+                )
 
             # Check if there is an attribute that defines the result dictionary keys
             if hasattr(self, "result_keys"):
@@ -442,8 +416,7 @@ def osipi_fit_full_volume(self, data, bvalues=None, **kwargs):
             for key in self.result_keys:
                 results[key] = np.empty(list(data.shape[:-1]))
             # no normalisation as volume algorithms may not want normalized signals...
-            args = [data, use_bvalues]
-            fit = self.ivim_fit_full_volume(*args, **kwargs) # Assume this is a dict with an array per key representing the parametric maps
+            fit = self.ivim_fit_full_volume(data, **kwargs) # Assume this is a dict with an array per key representing the parametric maps
             for key in list(fit.keys()):
                 results[key] = fit[key]
 
@@ -572,9 +545,9 @@ def osipi_author(self):
         """Author identification"""
         return ''
     
-    def osipi_simple_bias_and_RMSE_test(self, SNR, bvalues, f, Dstar, D, noise_realizations=100):
-        # Generate signal
-        bvalues = np.asarray(bvalues)
+    def osipi_simple_bias_and_RMSE_test(self, SNR, f, Dstar, D, noise_realizations=100):
+        # Generate signal using b-values from initialization
+        bvalues = np.asarray(self.bvalues)
         signals = f*np.exp(-bvalues*Dstar) + (1-f)*np.exp(-bvalues*D)
         
         f_estimates = np.zeros(noise_realizations)
@@ -587,7 +560,7 @@ def osipi_simple_bias_and_RMSE_test(self, SNR, bvalues, f, Dstar, D, noise_reali
             
             # Perform fit with the noised signal
             # f_estimates[i], Dstar_estimates[i], D_estimates[i] = self.D_and_Ds_swap(self.ivim_fit(noised_signal, bvalues))
-            result = self.ivim_fit(noised_signal, bvalues)
+            result = self.ivim_fit(noised_signal)
             f_estimates[i], Dstar_estimates[i], D_estimates[i] = result['f'], result['Dp'], result['D']
             
         # Calculate bias
diff --git a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py
index 6df9cbbb..3debf519 100644
--- a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py
+++ b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py
@@ -29,14 +29,14 @@ def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
 
     #model = ETP_SRI_LinearFitting(thresholds=[200])
     if kwargs:
-        results = model.osipi_fit(signals, bvalues, **kwargs)
+        results = model.osipi_fit(signals, **kwargs)
     else:
-        results = model.osipi_fit(signals, bvalues)
+        results = model.osipi_fit(signals)
     print(results)
     #test = model.osipi_simple_bias_and_RMSE_test(SNR=20, bvalues=bvalues, f=0.1, Dstar=0.03, D=0.001, noise_realizations=10)
     
 #model1 = ETP_SRI_LinearFitting(thresholds=[200])
-model2 = IAR_LU_biexp()
+model2 = IAR_LU_biexp(bvalues=np.array([0, 50, 0, 100, 0, 150, 200, 500, 0, 800]))
 #model2 = IAR_LU_modified_mix()
 #model2 = OGC_AmsterdamUMC_biexp()
 
diff --git a/tests/IVIMmodels/unit_tests/test_ivim_fit.py b/tests/IVIMmodels/unit_tests/test_ivim_fit.py
index 9d1b74d9..5e7e61a6 100644
--- a/tests/IVIMmodels/unit_tests/test_ivim_fit.py
+++ b/tests/IVIMmodels/unit_tests/test_ivim_fit.py
@@ -51,9 +51,9 @@ def test_ivim_fit_saved(data_ivim_fit_saved, eng, request, record_property):
         test_bounds = {"S0" : [0.7, 1.3], "f" : [0, 1.0], "Dp" : [0.01, 0.2], "D" : [0, 0.005]}
     else:
         test_bounds = {"S0" : [0.7, 1.3], "f" : [0, 1.0], "Dp" : [0.005, 0.2], "D" : [0, 0.005]}
-    fit = OsipiBase(algorithm=algorithm, bounds=test_bounds, **kwargs)
+    fit = OsipiBase(algorithm=algorithm, bvalues=bvals, bounds=test_bounds, **kwargs)
     start_time = time.time()  # Record the start time
-    fit_result = fit.osipi_fit(signal, bvals)
+    fit_result = fit.osipi_fit(signal)
     elapsed_time = time.time() - start_time  # Calculate elapsed time
     def to_list_if_needed(value):
         return value.tolist() if isinstance(value, np.ndarray) else value
@@ -123,10 +123,10 @@ def test_bounds(bound_input, eng, request):
     #bounds = ([0.0008, 0.2, 0.01, 1.1], [0.0012, 0.3, 0.02, 1.3])
     bounds = {"S0" : [1.1, 1.3], "f" : [0.2, 0.3], "Dp" : [0.01, 0.02], "D" : [0.0008, 0.0012]}
     # deliberately have silly bounds to see whether they are used
-    fit = OsipiBase(algorithm=algorithm, bounds=bounds, initial_guess={"S0" : 1.2, "f" : 0.25, "Dp" : 0.015, "D" : 0.001}, **kwargs)
+    fit = OsipiBase(algorithm=algorithm, bvalues=bvals, bounds=bounds, initial_guess={"S0" : 1.2, "f" : 0.25, "Dp" : 0.015, "D" : 0.001}, **kwargs)
     if fit.use_bounds["f"] or fit.use_bounds["D"] or fit.use_bounds["Dp"]:
         signal = signal_helper(data["data"])
-        fit_result = fit.osipi_fit(signal, bvals)
+        fit_result = fit.osipi_fit(signal)
         eps=1e-10 # without this margin it can cause floating point failures on mac systems
         if fit.use_bounds["D"]:
             assert bounds["D"][0]-eps <= fit_result['D'] <= bounds["D"][1]+eps,  f"Result {fit_result['D']} out of bounds for data: {name}"
@@ -185,10 +185,10 @@ def test_volume(algorithmlist,eng, threeddata):
         kwargs = {'bvalues': bvals}
     else:
         kwargs={}
-    fit = OsipiBase(algorithm=algorithm,**kwargs)
+    fit = OsipiBase(algorithm=algorithm, bvalues=bvals, **kwargs)
     if hasattr(fit, 'ivim_fit_full_volume') and callable(getattr(fit, 'ivim_fit_full_volume')):
         start_time = time.time()  # Record the start time
-        fit_result = fit.osipi_fit_full_volume(data, bvals)
+        fit_result = fit.osipi_fit_full_volume(data)
         elapsed_time2 = time.time() - start_time  # Calculate elapsed time
         print('time elaapsed is '+str(elapsed_time2))
         assert np.shape(fit_result['D'])[0] == np.shape(data)[0]
@@ -218,17 +218,17 @@ def test_parallel(algorithmlist,eng,threeddata):
     data[invalid_mask,:] = np.nan
     print('testing ' + str(np.sum(~invalid_mask)) + ' voxels of a matrix size ' + str(np.shape(data)))
 
-    fit = OsipiBase(algorithm=algorithm,**kwargs)
+    fit = OsipiBase(algorithm=algorithm, bvalues=bvals, **kwargs)
 
     def dummy_task(x):
         return x
 
     start_time = time.time()  # Record the start time
-    fit_result = fit.osipi_fit(data, bvals,njobs=1)
+    fit_result = fit.osipi_fit(data, njobs=1)
     time_serial = time.time() - start_time  # Calculate elapsed time
     Parallel(n_jobs=2)(delayed(dummy_task)(i) for i in range(8)) #github actions only supports 2 cores
     start_time = time.time()  # Record the start time
-    fit_result2 = fit.osipi_fit(data, bvals,njobs=2)
+    fit_result2 = fit.osipi_fit(data, njobs=2)
     time_parallel= time.time() - start_time  # Calculate elapsed time
 
     print('singular took '+str(time_serial)+' seconds, and parallel took '+str(time_parallel)+' seconds')
@@ -259,7 +259,7 @@ def test_deep_learning_algorithms(deep_learning_algorithms, record_property):
     fit = OsipiBase(bvalues=bvals, algorithm=algorithm, bounds={"S0" : [0, 2], "f" : [0, 1], "Dp" : [0.005, 0.2], "D" : [0, 0.005]}, **kwargs)
 
     array_2d = np.array([dat["data"] for _, dat in data.items()])
-    fit_result = fit.osipi_fit_full_volume(array_2d, bvals)
+    fit_result = fit.osipi_fit_full_volume(array_2d)
 
     errors = []  # Collect all assertion errors
 
diff --git a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py
index 18d176a1..e0ede7da 100644
--- a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py
+++ b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py
@@ -26,7 +26,7 @@ def test_generated(algorithmlist, ivim_data, SNR, rtol, atol, fit_count, rician_
     D = data["D"]
     f = data["f"]
     Dp = data["Dp"]
-    fit = OsipiBase(algorithm=ivim_algorithm)
+    fit = OsipiBase(algorithm=ivim_algorithm, bvalues=bvals)
     # here is a prior
     if use_prior and hasattr(fit, "supported_priors") and fit.supported_priors:
         prior = [rng.normal(D, D/3, 10), rng.normal(f, f/3, 10), rng.normal(Dp, Dp/3, 10), rng.normal(1, 1/3, 10)]
@@ -39,7 +39,7 @@ def test_generated(algorithmlist, ivim_data, SNR, rtol, atol, fit_count, rician_
         # else:
         #     signal = data["data"]
         start_time = datetime.datetime.now()
-        fit_result = fit.osipi_fit(signal, bvals) #, prior_in=prior
+        fit_result = fit.osipi_fit(signal) #, prior_in=prior
         time_delta += datetime.datetime.now() - start_time
         if save_file is not None:
             save_file.writerow([ivim_algorithm, name, SNR, idx, f, Dp, D, fit_result["f"], fit_result["Dp"], fit_result["D"], *signal])