diff --git a/main/.buildinfo b/main/.buildinfo
index 4a90e3e8..5a55549a 100644
--- a/main/.buildinfo
+++ b/main/.buildinfo
@@ -1,4 +1,4 @@
# Sphinx build info version 1
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: 93cb2caf85b8568b4f9f9e405da78ea3
+config: ef67081b51e860f83bae226709f7d9ec
tags: 645f666f9bcd5a90fca523b33c5a78b7
diff --git a/main/_modules/eddymotion/cli/parser.html b/main/_modules/eddymotion/cli/parser.html
index bb0f2601..a1baa916 100644
--- a/main/_modules/eddymotion/cli/parser.html
+++ b/main/_modules/eddymotion/cli/parser.html
@@ -5,7 +5,7 @@
- eddymotion.cli.parser - eddymotion 24.0.0.dev54 documentation
+ eddymotion.cli.parser - eddymotion 24.0.0.dev61 documentation
@@ -179,7 +179,7 @@
fromeddymotion.exceptionsimportModelNotFittedError
+DEFAULT_MIN_S0=1e-5
+"""Minimum value when considering the :math:`S_{0}` DWI signal."""
+
+DEFAULT_MAX_S0=1.0
+"""Maximum value when considering the :math:`S_{0}` DWI signal."""
+
+DEFAULT_MAX_BVALUE=1000
+"""Maximum allowed value for the b-value."""
+
+DEFAULT_LOWB_THRESHOLD=50
+"""The lower bound for the b-value so that the orientation is considered a DW volume."""
+
+DEFAULT_HIGHB_THRESHOLD=10000
+"""A b-value cap for DWI data."""
+
+DEFAULT_CLIP_PERCENTILE=75
+"""Upper percentile threshold for intensity clipping."""
+
+DEFAULT_TIMEFRAME_MIDPOINT_TOL=1e-2
+"""Time frame tolerance in seconds."""
+
def_exec_fit(model,data,chunk=None):retval=model.fit(data)returnretval,chunk
-def_exec_predict(model,gradient,chunk=None,**kwargs):
+def_exec_predict(model,chunk=None,**kwargs):"""Propagate model parameters and call predict."""
- returnnp.squeeze(model.predict(gradient,S0=kwargs.pop("S0",None))),chunk
+ returnnp.squeeze(model.predict(**kwargs)),chunk
+[docs]
+ deffit(self,data,**kwargs):
+"""Abstract member signature of fit()."""
+ raiseNotImplementedError("Cannot call fit() on a BaseModel instance.")
+
+
+
+[docs]
+ defpredict(self,*args,**kwargs):
+"""Abstract member signature of predict()."""
+ raiseNotImplementedError("Cannot call predict() on a BaseModel instance.")
+
+
+
+
+
+[docs]
+classPETModel(BaseModel):
+"""A PET imaging realignment model based on B-Spline approximation."""
+
+ __slots__=("_t","_x","_xlim","_order","_coeff","_n_ctrl")
+
+ def__init__(self,timepoints=None,xlim=None,n_ctrl=None,order=3,**kwargs):
+"""
+ Create the B-Spline interpolating matrix.
+
+ Parameters:
+ -----------
+ timepoints : :obj:`list`
+ The timing (in sec) of each PET volume.
+ E.g., ``[15., 45., 75., 105., 135., 165., 210., 270., 330.,
+ 420., 540., 750., 1050., 1350., 1650., 1950., 2250., 2550.]``
+
+ n_ctrl : :obj:`int`
+ Number of B-Spline control points. If `None`, then one control point every
+ six timepoints will be used. The less control points, the smoother is the
+ model.
+
+ """
+ super.__init__(**kwargs)
+
+ iftimepointsisNoneorxlimisNone:
+ raiseTypeError("timepoints must be provided in initialization")
+
+ self._order=order
+
+ self._x=np.array(timepoints,dtype="float32")
+ self._xlim=xlim
+
+ ifself._x[0]<DEFAULT_TIMEFRAME_MIDPOINT_TOL:
+ raiseValueError("First frame midpoint should not be zero or negative")
+ ifself._x[-1]>(self._xlim-DEFAULT_TIMEFRAME_MIDPOINT_TOL):
+ raiseValueError("Last frame midpoint should not be equal or greater than duration")
+
+ # Calculate index coordinates in the B-Spline grid
+ self._n_ctrl=n_ctrlor(len(timepoints)//4)+1
+
+ # B-Spline knots
+ self._t=np.arange(-3,float(self._n_ctrl)+4,dtype="float32")
+
+ self._coeff=None
+
+ @property
+ defis_fitted(self):
+ returnself._coeffisnotNone
+
+
+[docs]
+ deffit(self,data,**kwargs):
+"""Fit the model."""
+ fromscipy.interpolateimportBSpline
+ fromscipy.sparse.linalgimportcg
+
+ n_jobs=kwargs.pop("n_jobs",None)or1
+
+ timepoints=kwargs.get("timepoints",None)orself._x
+ x=(np.array(timepoints,dtype="float32")/self._xlim)*self._n_ctrl
+
+ self._datashape=data.shape[:3]
+
+ # Convert data into V (voxels) x T (timepoints)
+ data=data.reshape((-1,data.shape[-1]))ifself._maskisNoneelsedata[self._mask]
+
+ # A.shape = (T, K - 4); T= n. timepoints, K= n. knots (with padding)
+ A=BSpline.design_matrix(x,self._t,k=self._order)
+ AT=A.T
+ ATdotA=AT@A
+
+ # One single CPU - linear execution (full model)
+ ifn_jobs==1:
+ self._coeff=np.array([cg(ATdotA,AT@v)[0]forvindata])
+ return
+
+ # Parallelize process with joblib
+ withParallel(n_jobs=n_jobs)asexecutor:
+ results=executor(delayed(cg)(ATdotA,AT@v)forvindata)
+
+ self._coeff=np.array([r[0]forrinresults])
+
+
+
+[docs]
+ defpredict(self,index=None,**kwargs):
+"""Return the corrected volume using B-spline interpolation."""
+ fromscipy.interpolateimportBSpline
+
+ ifindexisNone:
+ raiseValueError("A timepoint index to be simulated must be provided.")
+
+ ifnotself._is_fitted:
+ raiseModelNotFittedError(f"{type(self).__name__} must be fitted before predicting")
+
+ # Project sample timing into B-Spline coordinates
+ x=(index/self._xlim)*self._n_ctrl
+ A=BSpline.design_matrix(x,self._t,k=self._order)
+
+ # A is 1 (num. timepoints) x C (num. coeff)
+ # self._coeff is V (num. voxels) x K - 4
+ predicted=np.squeeze(A@self._coeff.T)
+
+ ifself._maskisNone:
+ returnpredicted.reshape(self._datashape)
+
+ retval=np.zeros(self._datashape,dtype="float32")
+ retval[self._mask]=predicted
+ returnretval
+
+
+
+
+
+[docs]
+classBaseDWIModel(BaseModel):
+"""Interface and default methods for DWI models."""
+
+ __slots__=(
+ "_gtab",
+ "_S0",
+ "_b_max",
+ "_model_class",# Defining a model class, DIPY models are instantiated automagically
+ "_modelargs",
+ )
+
+ def__init__(self,gtab,S0=None,b_max=None,**kwargs):
+"""Initialization.
+
+ Parameters
+ ----------
+ gtab : :obj:`numpy.ndarray`
+ An :math:`N \times 4` table, where rows (*N*) are diffusion gradients and
+ columns are b-vector components and corresponding b-value, respectively.
+ S0 : :obj:`numpy.ndarray`
+ :math:`S_{0}` signal.
+ b_max : :obj:`int`
+ Maximum value to cap b-values.
+
+ """
+
+ super().__init__(**kwargs)
+
# Setup B0 mapself._S0=NoneifS0isnotNone:self._S0=np.clip(S0.astype("float32")/S0.max(),
- a_min=1e-5,
- a_max=1.0,
+ a_min=DEFAULT_MIN_S0,
+ a_max=DEFAULT_MAX_S0,)
- # Setup brain mask
- self._mask=mask
- ifmaskisNoneandS0isnotNone:
- self._mask=self._S0>np.percentile(self._S0,35)
-
# Cap b-values, if requested
+ self._gtab=gtabself._b_max=None
- ifb_maxandb_max>1000:
+ ifb_maxandb_max>DEFAULT_MAX_BVALUE:# Saturate b-values at b_max, since signal stops dropping
- gtab[-1,gtab[-1]>b_max]=b_max
+ self._gtab[-1,self._gtab[-1]>b_max]=b_max# A possibly good alternative is completely remove very high b-values# bval_mask = gtab[-1] < b_max# data = data[..., bval_mask]
@@ -429,25 +620,19 @@
+[docs]
+ defpredict(self,gradient=None,**kwargs):"""Predict asynchronously chunk-by-chunk the diffusion signal."""
+ ifgradientisNone:
+ raiseValueError("A gradient to be simulated (b-vector, b-value) must be provided")
+
ifnotself._is_fitted:raiseModelNotFittedError(f"{type(self).__name__} must be fitted before predicting")
- ifself._b_maxisnotNone:
- gradient[-1]=min(gradient[-1],self._b_max)
+ gradient=np.array(gradient)# Tuples are unmutable
+
+ # Cap the b-value if b_max is defined
+ gradient[-1]=min(gradient[-1],self._b_maxorgradient[-1])gradient=_rasb2dipy(gradient)
@@ -505,7 +695,7 @@
[docs]
- defpredict(self,gradient,**kwargs):
+ defpredict(self,*_,**kwargs):"""Return the *b=0* map."""# No need to check fit (if not fitted, has raised already)
@@ -570,7 +762,7 @@
Source code for eddymotion.model.base
[docs]
-classAverageDWModel:
+classAverageDWModel(BaseDWIModel):"""A trivial model that returns an average map."""__slots__=("_data","_th_low","_th_high","_bias","_stat","_is_fitted")
@@ -581,24 +773,24 @@
Source code for eddymotion.model.base
Parameters ----------
- gtab : :obj:`~numpy.ndarray`
- An :math:`N \times 4` table, where rows (*N*) are diffusion gradients and
- columns are b-vector components and corresponding b-value, respectively.
- th_low : :obj:`~numbers.Number`
+ th_low : :obj:`numbers.Number` A lower bound for the b-value corresponding to the diffusion weighted images that will be averaged.
- th_high : :obj:`~numbers.Number`
+ th_high : :obj:`numbers.Number` An upper bound for the b-value corresponding to the diffusion weighted images that will be averaged. bias : :obj:`bool` Whether the overall distribution of each diffusion weighted image will be
- standardized and centered around the global 75th percentile.
+ standardized and centered around the
+ :data:`src.eddymotion.model.base.DEFAULT_CLIP_PERCENTILE` percentile. stat : :obj:`str` Whether the summary statistic to apply is ``"mean"`` or ``"median"``. """
- self._th_low=kwargs.get("th_low",50)
- self._th_high=kwargs.get("th_high",10000)
+ super().__init__(**kwargs)
+
+ self._th_low=kwargs.get("th_low",DEFAULT_LOWB_THRESHOLD)
+ self._th_high=kwargs.get("th_high",DEFAULT_HIGHB_THRESHOLD)self._bias=kwargs.get("bias",True)self._stat=kwargs.get("stat","median")self._data=None
@@ -608,7 +800,10 @@
Source code for eddymotion.model.base
[docs]deffit(self,data,**kwargs):"""Calculate the average."""
- gtab=kwargs.pop("gtab",None)
+
+ if(gtab:=kwargs.pop("gtab",None))isNone:
+ raiseValueError("A gradient table must be provided.")
+
# Select the interval of b-values for which DWIs will be averagedb_mask=(((gtab[3]>=self._th_low)&(gtab[3]<=self._th_high))
@@ -620,7 +815,7 @@
Source code for eddymotion.model.base
# Regress out global signal differencesifself._bias:centers=np.median(shells,axis=(0,1,2))
- reference=np.percentile(centers[centers>=1.0],75)
+ reference=np.percentile(centers[centers>=1.0],DEFAULT_CLIP_PERCENTILE)centers[centers<1.0]=referencedrift=reference/centersshells=shells*drift
@@ -629,8 +824,7 @@
Source code for eddymotion.model.base
avg_func=np.medianifself._stat=="median"elsenp.mean# Calculate the averageself._data=avg_func(shells,axis=-1)
-
- self._is_fitted=self._dataisnotNone
+ self._is_fitted=True
@property
@@ -639,7 +833,7 @@
Source code for eddymotion.model.base
[docs]
- defpredict(self,gradient,**kwargs):
+ defpredict(self,*_,**kwargs):"""Return the average map."""ifnotself._is_fitted:
@@ -650,121 +844,9 @@
Source code for eddymotion.model.base
-
-[docs]
-classPETModel:
-"""A PET imaging realignment model based on B-Spline approximation."""
-
- __slots__=("_t","_x","_xlim","_order","_coeff","_mask","_shape","_n_ctrl")
-
- def__init__(self,timepoints=None,xlim=None,n_ctrl=None,mask=None,order=3,**kwargs):
-"""
- Create the B-Spline interpolating matrix.
-
- Parameters:
- -----------
- timepoints : :obj:`list`
- The timing (in sec) of each PET volume.
- E.g., ``[15., 45., 75., 105., 135., 165., 210., 270., 330.,
- 420., 540., 750., 1050., 1350., 1650., 1950., 2250., 2550.]``
-
- n_ctrl : :obj:`int`
- Number of B-Spline control points. If `None`, then one control point every
- six timepoints will be used. The less control points, the smoother is the
- model.
-
- """
- iftimepointsisNoneorxlimisNone:
- raiseTypeError("timepoints must be provided in initialization")
-
- self._order=order
- self._mask=mask
-
- self._x=np.array(timepoints,dtype="float32")
- self._xlim=xlim
-
- ifself._x[0]<1e-2:
- raiseValueError("First frame midpoint should not be zero or negative")
- ifself._x[-1]>(self._xlim-1e-2):
- raiseValueError("Last frame midpoint should not be equal or greater than duration")
-
- # Calculate index coordinates in the B-Spline grid
- self._n_ctrl=n_ctrlor(len(timepoints)//4)+1
-
- # B-Spline knots
- self._t=np.arange(-3,float(self._n_ctrl)+4,dtype="float32")
-
- self._shape=None
- self._coeff=None
-
- @property
- defis_fitted(self):
- returnself._coeffisnotNone
-
-
-[docs]
- deffit(self,data,*args,**kwargs):
-"""Fit the model."""
- fromscipy.interpolateimportBSpline
- fromscipy.sparse.linalgimportcg
-
- n_jobs=kwargs.pop("n_jobs",None)or1
-
- timepoints=kwargs.get("timepoints",None)orself._x
- x=(np.array(timepoints,dtype="float32")/self._xlim)*self._n_ctrl
-
- self._shape=data.shape[:3]
-
- # Convert data into V (voxels) x T (timepoints)
- data=data.reshape((-1,data.shape[-1]))ifself._maskisNoneelsedata[self._mask]
-
- # A.shape = (T, K - 4); T= n. timepoints, K= n. knots (with padding)
- A=BSpline.design_matrix(x,self._t,k=self._order)
- AT=A.T
- ATdotA=AT@A
-
- # One single CPU - linear execution (full model)
- ifn_jobs==1:
- self._coeff=np.array([cg(ATdotA,AT@v)[0]forvindata])
- return
-
- # Parallelize process with joblib
- withParallel(n_jobs=n_jobs)asexecutor:
- results=executor(delayed(cg)(ATdotA,AT@v)forvindata)
-
- self._coeff=np.array([r[0]forrinresults])
-
-
-
-[docs]
- defpredict(self,timepoint,**kwargs):
-"""Return the *b=0* map."""
- fromscipy.interpolateimportBSpline
-
- ifnotself._is_fitted:
- raiseModelNotFittedError(f"{type(self).__name__} must be fitted before predicting")
-
- # Project sample timing into B-Spline coordinates
- x=(timepoint/self._xlim)*self._n_ctrl
- A=BSpline.design_matrix(x,self._t,k=self._order)
-
- # A is 1 (num. timepoints) x C (num. coeff)
- # self._coeff is V (num. voxels) x K - 4
- predicted=np.squeeze(A@self._coeff.T)
-
- ifself._maskisNone:
- returnpredicted.reshape(self._shape)
-
- retval=np.zeros(self._shape,dtype="float32")
- retval[self._mask]=predicted
- returnretval
-
-
-
-
[docs]
-classDTIModel(BaseModel):
+classDTIModel(BaseDWIModel):"""A wrapper of :obj:`dipy.reconst.dti.TensorModel`."""_modelargs=(
@@ -781,7 +863,7 @@
Source code for eddymotion.model.base
[docs]
-classDKIModel(BaseModel):
+classDKIModel(BaseDWIModel):"""A wrapper of :obj:`dipy.reconst.dki.DiffusionKurtosisModel`."""_modelargs=DTIModel._modelargs
@@ -837,7 +919,7 @@