diff --git a/requirements.txt b/requirements.txt index 29742e3..16244cb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ numpy<2 +super-ivim-dc==0.4.0 nibabel scipy torchio diff --git a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py index 455e639..d760982 100644 --- a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py +++ b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py @@ -5,7 +5,7 @@ class ASD_MemorialSloanKettering_QAMPER_IVIM(OsipiBase): """ - Bi-exponential fitting algorithm by Oliver Gurney-Champion, Amsterdam UMC + Bi-exponential fitting algorithm by Eve LoCastro and Amita Shukla-Dave, Memorial Sloan Kettering """ # I'm thinking that we define default attributes for each submission like this diff --git a/src/standardized/TCML_TechnionIIT_SLS.py b/src/standardized/TCML_TechnionIIT_SLS.py new file mode 100644 index 0000000..9021a1d --- /dev/null +++ b/src/standardized/TCML_TechnionIIT_SLS.py @@ -0,0 +1,95 @@ +from src.wrappers.OsipiBase import OsipiBase +from super_ivim_dc.source.Classsic_ivim_fit import IVIM_fit_sls +import numpy as np + +class TCML_TechnionIIT_SLS(OsipiBase): + """ + TCML_TechnionIIT_lsqBOBYQA fitting algorithm by Moti Freiman and Noam Korngut, TechnionIIT + """ + + # I'm thinking that we define default attributes for each submission like this + # And in __init__, we can call the OsipiBase control functions to check whether + # the user inputs fulfil the requirements + + # Some basic stuff that identifies the algorithm + id_author = "Moti Freiman and Noam Korngut, TechnIIT" + id_algorithm_type = "Bi-exponential fit, segmented fitting" + id_return_parameters = "f, D*, D, S0" + id_units = "seconds per milli metre squared or milliseconds per micro metre squared" + + # Algorithm requirements + required_bvalues = 4 + required_thresholds = [0,0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds + required_bounds = False + required_bounds_optional = True # Bounds may not be required but are optional + required_initial_guess = False + required_initial_guess_optional = True + accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most? + + + # Supported inputs in the standardized class + supported_bounds = True + supported_initial_guess = True + supported_thresholds = True + + def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True): + """ + Everything this algorithm requires should be implemented here. + Number of segmentation thresholds, bounds, etc. + + Our OsipiBase object could contain functions that compare the inputs with + the requirements. + """ + super(TCML_TechnionIIT_SLS, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess) + self.fit_least_squares = IVIM_fit_sls + self.fitS0=fitS0 + self.initialize(bounds, initial_guess, fitS0,thresholds) + + def initialize(self, bounds, initial_guess, fitS0,thresholds): + if bounds is None: + print('warning, no bounds were defined, so default bounds are used of ([0.0003, 0.001, 0.009, 0],[0.008, 0.5,0.04, 3])') + self.bounds = ([0.0003, 0.001, 0.009, 0],[0.008, 0.5,0.04, 3]) + else: + bounds=bounds + self.bounds = bounds + if initial_guess is None: + print('warning, no initial guesses were defined, so default bounds are used of [0.001, 0.1, 0.01, 1]') + self.initial_guess = [0.001, 0.1, 0.01, 1] # D, Dp, f, S0 + else: + self.initial_guess = initial_guess + self.use_initial_guess = True + self.fitS0=fitS0 + self.use_initial_guess = True + self.use_bounds = True + if thresholds is None: + self.thresholds = 150 + print('warning, no thresholds were defined, so default bounds are used of 150') + else: + self.thresholds = thresholds + + def ivim_fit(self, signals, bvalues, **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) + bounds=np.array(self.bounds) + bounds=[bounds[0][[0, 2, 1, 3]], bounds[1][[0, 2, 1, 3]]] + initial_guess = np.array(self.initial_guess) + initial_guess = initial_guess[[0, 2, 1, 3]] + + + fit_results = self.fit_least_squares(1 ,np.array(signals)[:,np.newaxis],bvalues, bounds,initial_guess,self.thresholds) + + results = {} + results["D"] = fit_results[0] + results["f"] = fit_results[2] + results["Dp"] = fit_results[1] + + return results \ No newline at end of file diff --git a/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py b/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py new file mode 100644 index 0000000..169f8bd --- /dev/null +++ b/src/standardized/TCML_TechnionIIT_lsqBOBYQA.py @@ -0,0 +1,90 @@ +from src.wrappers.OsipiBase import OsipiBase +from super_ivim_dc.source.Classsic_ivim_fit import fit_least_squers_BOBYQA +import numpy as np + +class TCML_TechnionIIT_lsqBOBYQA(OsipiBase): + """ + TCML_TechnionIIT_lsqBOBYQA fitting algorithm by Moti Freiman and Noam Korngut, TechnionIIT + """ + + # I'm thinking that we define default attributes for each submission like this + # And in __init__, we can call the OsipiBase control functions to check whether + # the user inputs fulfil the requirements + + # Some basic stuff that identifies the algorithm + id_author = "Moti Freiman and Noam Korngut, TechnIIT" + id_algorithm_type = "Bi-exponential fit, BOBYQO" + id_return_parameters = "f, D*, D, S0" + id_units = "seconds per milli metre squared or milliseconds per micro metre squared" + + # Algorithm requirements + required_bvalues = 4 + required_thresholds = [0, + 0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds + required_bounds = False + required_bounds_optional = True # Bounds may not be required but are optional + required_initial_guess = False + required_initial_guess_optional = True + accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most? + + + # Supported inputs in the standardized class + supported_bounds = True + supported_initial_guess = True + supported_thresholds = False + + def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True): + """ + Everything this algorithm requires should be implemented here. + Number of segmentation thresholds, bounds, etc. + + Our OsipiBase object could contain functions that compare the inputs with + the requirements. + """ + super(TCML_TechnionIIT_lsqBOBYQA, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess) + self.fit_least_squares = fit_least_squers_BOBYQA + self.fitS0=fitS0 + self.initialize(bounds, initial_guess, fitS0) + + def initialize(self, bounds, initial_guess, fitS0): + if bounds is None: + print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]') + self.bounds = ([0.0003, 0.001, 0.009, 0],[0.01, 0.5,0.04, 3]) + else: + bounds=bounds + self.bounds = bounds + if initial_guess is None: + print('warning, no initial guesses were defined, so default bounds are used of [0.001, 0.001, 0.01, 1]') + self.initial_guess = [0.001, 0.1, 0.01, 1] # D, Dp, f, S0 + else: + self.initial_guess = initial_guess + self.use_initial_guess = True + self.fitS0=fitS0 + self.use_initial_guess = True + self.use_bounds = True + + def ivim_fit(self, signals, bvalues, **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) + bounds=np.array(self.bounds) + bounds=[bounds[0][[0, 2, 1, 3]], bounds[1][[0, 2, 1, 3]]] + initial_guess = np.array(self.initial_guess) + initial_guess = initial_guess[[0, 2, 1, 3]] + + fit_results = self.fit_least_squares(1 ,bvalues, np.array(signals)[:,np.newaxis], bounds, initial_guess) + + results = {} + results["D"] = fit_results[0] + results["f"] = fit_results[2] + results["Dp"] = fit_results[1] + + return results \ No newline at end of file diff --git a/src/standardized/TCML_TechnionIIT_lsqlm.py b/src/standardized/TCML_TechnionIIT_lsqlm.py new file mode 100644 index 0000000..86af016 --- /dev/null +++ b/src/standardized/TCML_TechnionIIT_lsqlm.py @@ -0,0 +1,81 @@ +from src.wrappers.OsipiBase import OsipiBase +from super_ivim_dc.source.Classsic_ivim_fit import fit_least_squares_lm +import numpy as np + +class TCML_TechnionIIT_lsqlm(OsipiBase): + """ + TCML_TechnionIIT_lsqlm fitting algorithm by Moti Freiman and Noam Korngut, TechnionIIT + """ + + # I'm thinking that we define default attributes for each submission like this + # And in __init__, we can call the OsipiBase control functions to check whether + # the user inputs fulfil the requirements + + # Some basic stuff that identifies the algorithm + id_author = "Moti Freiman and Noam Korngut, TechnIIT" + id_algorithm_type = "Bi-exponential fit" + id_return_parameters = "f, D*, D, S0" + id_units = "seconds per milli metre squared or milliseconds per micro metre squared" + + # Algorithm requirements + required_bvalues = 4 + required_thresholds = [0, + 0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds + required_bounds = False + required_bounds_optional = False # Bounds may not be required but are optional + required_initial_guess = False + required_initial_guess_optional = True + accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most? + + + # Supported inputs in the standardized class + supported_bounds = False + supported_initial_guess = True + supported_thresholds = False + + def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True): + """ + Everything this algorithm requires should be implemented here. + Number of segmentation thresholds, bounds, etc. + + Our OsipiBase object could contain functions that compare the inputs with + the requirements. + """ + super(TCML_TechnionIIT_lsqlm, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess) + self.fit_least_squares = fit_least_squares_lm + self.fitS0=fitS0 + self.initialize(bounds, initial_guess, fitS0) + + def initialize(self, bounds, initial_guess, fitS0): + if initial_guess is None: + print('warning, no initial guesses were defined, so default bounds are used of [0.001, 0.1, 0.01, 1]') + self.initial_guess = [0.001, 0.1, 0.01, 1] + else: + self.initial_guess = initial_guess + self.use_initial_guess = True + self.fitS0=fitS0 + self.use_initial_guess = True + self.use_bounds = False + + def ivim_fit(self, signals, bvalues, **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) + initial_guess = np.array(self.initial_guess) + initial_guess = initial_guess[[0, 2, 1, 3]] + fit_results = self.fit_least_squares(1 ,bvalues, np.array(signals)[:,np.newaxis], initial_guess) + + results = {} + results["D"] = fit_results[0] + results["f"] = fit_results[2] + results["Dp"] = fit_results[1] + + return results \ No newline at end of file diff --git a/src/standardized/TCML_TechnionIIT_lsqtrf.py b/src/standardized/TCML_TechnionIIT_lsqtrf.py new file mode 100644 index 0000000..bf7083d --- /dev/null +++ b/src/standardized/TCML_TechnionIIT_lsqtrf.py @@ -0,0 +1,89 @@ +from src.wrappers.OsipiBase import OsipiBase +from super_ivim_dc.source.Classsic_ivim_fit import fit_least_squares_trf +import numpy as np + +class TCML_TechnionIIT_lsqtrf(OsipiBase): + """ + TCML_TechnionIIT_lsqlm fitting algorithm by Moti Freiman and Noam Korngut, TechnionIIT + """ + + # I'm thinking that we define default attributes for each submission like this + # And in __init__, we can call the OsipiBase control functions to check whether + # the user inputs fulfil the requirements + + # Some basic stuff that identifies the algorithm + id_author = "Moti Freiman and Noam Korngut, TechnIIT" + id_algorithm_type = "Bi-exponential fit, Trust Region Reflective algorithm" + id_return_parameters = "f, D*, D, S0" + id_units = "seconds per milli metre squared or milliseconds per micro metre squared" + + # Algorithm requirements + required_bvalues = 4 + required_thresholds = [0, + 0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds + required_bounds = False + required_bounds_optional = True # Bounds may not be required but are optional + required_initial_guess = False + required_initial_guess_optional = True + accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most? + + + # Supported inputs in the standardized class + supported_bounds = True + supported_initial_guess = True + supported_thresholds = False + + def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True): + """ + Everything this algorithm requires should be implemented here. + Number of segmentation thresholds, bounds, etc. + + Our OsipiBase object could contain functions that compare the inputs with + the requirements. + """ + super(TCML_TechnionIIT_lsqtrf, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess) + self.fit_least_squares = fit_least_squares_trf + self.fitS0=fitS0 + self.initialize(bounds, initial_guess, fitS0) + + def initialize(self, bounds, initial_guess, fitS0): + if bounds is None: + print('warning, no bounds were defined, so default bounds are used of ([0.0003, 0.001, 0.009, 0],[0.008, 0.5,0.04, 3])') + self.bounds = ([0.0003, 0.001, 0.009, 0],[0.008, 0.5,0.04, 3]) + else: + bounds=bounds + self.bounds = bounds + if initial_guess is None: + print('warning, no initial guesses were defined, so default bounds are used of [0.001, 0.1, 0.01, 1]') + self.initial_guess = [0.001, 0.1, 0.01, 1] # D, Dp, f, S0 + else: + self.initial_guess = initial_guess + self.use_initial_guess = True + self.fitS0=fitS0 + self.use_initial_guess = True + self.use_bounds = True + + def ivim_fit(self, signals, bvalues, **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) + bounds=np.array(self.bounds) + bounds=[bounds[0][[0, 2, 1, 3]], bounds[1][[0, 2, 1, 3]]] + initial_guess = np.array(self.initial_guess) + initial_guess = initial_guess[[0, 2, 1, 3]] + fit_results = self.fit_least_squares(1 ,bvalues, np.array(signals)[:,np.newaxis], bounds,initial_guess) + + results = {} + results["D"] = fit_results[0] + results["f"] = fit_results[2] + results["Dp"] = fit_results[1] + + return results \ No newline at end of file diff --git a/src/wrappers/OsipiBase.py b/src/wrappers/OsipiBase.py index 6c5d380..2875ab9 100644 --- a/src/wrappers/OsipiBase.py +++ b/src/wrappers/OsipiBase.py @@ -106,7 +106,7 @@ def osipi_fit(self, data, bvalues=None, **kwargs): for ijk in tqdm(np.ndindex(data.shape[:-1]), total=np.prod(data.shape[:-1])): args = [data[ijk], use_bvalues] - fit = self.ivim_fit(*args, **kwargs) # For single voxel fits, we assume this is a dict with a float value per key. + fit = self.D_and_Ds_swap(self.ivim_fit(*args, **kwargs)) # For single voxel fits, we assume this is a dict with a float value per key. for key in list(fit.keys()): results[key][ijk] = fit[key] @@ -282,7 +282,7 @@ def osipi_simple_bias_and_RMSE_test(self, SNR, bvalues, f, Dstar, D, noise_reali noised_signal = np.array([norm.rvs(signal, sigma) for signal in signals]) # Perform fit with the noised signal - f_estimates[i], Dstar_estimates[i], D_estimates[i] = self.ivim_fit(noised_signal, bvalues) + f_estimates[i], Dstar_estimates[i], D_estimates[i] = self.D_and_Ds_swap(self.ivim_fit(noised_signal, bvalues)) # Calculate bias f_bias = np.mean(f_estimates) - f @@ -299,3 +299,10 @@ def osipi_simple_bias_and_RMSE_test(self, SNR, bvalues, f, Dstar, D, noise_reali print(f"D bias:\t{D_bias}\nD RMSE:\t{D_RMSE}") + def D_and_Ds_swap(self,results): + if results['D']>results['Dp']: + D=results['Dp'] + results['Dp']=results['D'] + results['D']=D + results['f']=1-results['f'] + return results \ No newline at end of file diff --git a/tests/IVIMmodels/unit_tests/algorithms.json b/tests/IVIMmodels/unit_tests/algorithms.json index fa90ce5..2362330 100644 --- a/tests/IVIMmodels/unit_tests/algorithms.json +++ b/tests/IVIMmodels/unit_tests/algorithms.json @@ -1,5 +1,7 @@ { "algorithms": [ + "TCML_TechnionIIT_lsqlm", + "TCML_TechnionIIT_lsqtrf", "ASD_MemorialSloanKettering_QAMPER_IVIM", "ETP_SRI_LinearFitting", "IAR_LU_biexp", diff --git a/tests/IVIMmodels/unit_tests/test_ivim_fit.py b/tests/IVIMmodels/unit_tests/test_ivim_fit.py index 9ed85b7..a38d1ba 100644 --- a/tests/IVIMmodels/unit_tests/test_ivim_fit.py +++ b/tests/IVIMmodels/unit_tests/test_ivim_fit.py @@ -42,7 +42,9 @@ def test_ivim_fit_saved(data_ivim_fit_saved, eng, request, record_property): request.node.add_marker(mark) signal = signal_helper(data["data"]) tolerances = tolerances_helper(tolerances, data) - fit = OsipiBase(algorithm=algorithm, **kwargs) + fit = OsipiBase(algorithm=algorithm, **kwargs) + if fit.use_bounds: + fit.bounds = ([0, 0, 0.005, 0.7], [0.005, 1.0, 0.2, 1.3]) start_time = time.time() # Record the start time fit_result = fit.osipi_fit(signal, bvals) elapsed_time = time.time() - start_time # Calculate elapsed time @@ -61,6 +63,8 @@ def to_list_if_needed(value): "atol": tolerances["atol"] } record_property('test_data', test_result) + if (data['f'] > 0.99 or data['f'] < 0.01) and not fit.use_bounds: #in these cases there are multiple solutions (D can become D*, f will be 0 and D* can be anything. This will be a good description of the data, so technically not a fail + return npt.assert_allclose(fit_result['f'],data['f'], rtol=tolerances["rtol"]["f"], atol=tolerances["atol"]["f"]) if data['f']<0.80: # we need some signal for D to be detected npt.assert_allclose(fit_result['D'],data['D'], rtol=tolerances["rtol"]["D"], atol=tolerances["atol"]["D"])