diff --git a/.gitignore b/.gitignore index 7e7c4a57..08036081 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,4 @@ __pycache__/ *.cpp dist/ -trk2dictionary.c \ No newline at end of file +trk2dictionary.c diff --git a/CHANGELOG.md b/CHANGELOG.md index 4224780e..269c22c1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,16 @@ # Change Log All notable changes to COMMIT will be documented in this file. +## [1.5.0] - 2021-02-16 + +### Added +- core.pyx: Add to the function build_operator the parameter tikhonov_equalizer and deriv_matrix +- operator.pyx: Extend Ax and A'y multiplications when Tikhonov regularization is enabled +- operator_withLUC.c: Add C functions to extend Ax and A'y multiplications when Tikhonov regularization is enabled + +### Changed +- core.pyx: The function get_y returns a larger vector y when Tikhonov regularization is enabled + ## [1.4.5] - 2021-02-08 ### Fixed diff --git a/commit/core.pyx b/commit/core.pyx index dba0be2f..511eb257 100755 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -639,7 +639,7 @@ cdef class Evaluation : LOG( ' [ %.1f seconds ]' % ( time.time() - tic ) ) - def build_operator( self, build_dir=None ) : + def build_operator( self, build_dir=None, tikhonov_lambda=0, tikhonov_matrix=None ) : """Compile/build the operator for computing the matrix-vector multiplications by A and A' using the informations from self.DICTIONARY, self.KERNELS and self.THREADS. NB: needs to call this function to update pointers to data structures in case @@ -652,6 +652,11 @@ cdef class Evaluation : If None (default), they will end up in the .pyxbld directory in the user’s home directory. If using this option, it is recommended to use a temporary directory, quit your python console between each build, and delete the content of the temporary directory. + tikhonov_lambda: float + Tikhonov lambda + If a positive value is given, tikhonov_matrix must not be None + tikhonov_matrix: string + Tikhonov matrix """ if self.DICTIONARY is None : ERROR( 'Dictionary not loaded; call "load_dictionary()" first' ) @@ -659,6 +664,12 @@ cdef class Evaluation : ERROR( 'Response functions not generated; call "generate_kernels()" and "load_kernels()" first' ) if self.THREADS is None : ERROR( 'Threads not set; call "set_threads()" first' ) + if tikhonov_lambda < 0: + ERROR( 'Invalid lambda for Tikhonov regularization; value must be positive or zero' ) + if tikhonov_lambda > 0 and tikhonov_matrix == None: + ERROR( 'Tikhonov lambda given but Tikhonov matrix was not selected; add "tikhonov_matrix" parameter in "build_operator()"' ) + if tikhonov_lambda > 0 and tikhonov_matrix!='L1' and tikhonov_matrix!='L2' and tikhonov_matrix!='L1z' and tikhonov_matrix!='L2z': + ERROR( 'Invalid matrix selection for Tikhonov regularization term; check "tikhonov_matrix" parameter in "build_operator()"' ) if self.DICTIONARY['IC']['nF'] <= 0 : ERROR( 'No streamline found in the dictionary; check your data' ) @@ -710,7 +721,7 @@ cdef class Evaluation : else : reload( sys.modules['commit.operator.operator'] ) - self.A = sys.modules['commit.operator.operator'].LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS ) + self.A = sys.modules['commit.operator.operator'].LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, tikhonov_lambda, tikhonov_matrix ) LOG( ' [ %.1f seconds ]' % ( time.time() - tic ) ) @@ -724,7 +735,26 @@ cdef class Evaluation : ERROR( 'Dictionary not loaded; call "load_dictionary()" first' ) if self.niiDWI is None : ERROR( 'Data not loaded; call "load_data()" first' ) - return self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float64) + + y = self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float64) + + # extend y for the tikhonov regularization term + if self.A.tikhonov_lambda > 0: + if self.A.tikhonov_matrix == 'L1': + yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-1, dtype=np.float64) + elif self.A.tikhonov_matrix == 'L2': + yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-2, dtype=np.float64) + elif self.A.tikhonov_matrix == 'L1z': + yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]+1, dtype=np.float64) + elif self.A.tikhonov_matrix == 'L2z': + yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0] , dtype=np.float64) + else: + ERROR( 'Invalid matrix selection for Tikhonov regularization term; check "tikhonov_matrix" parameter in "build_operator()"' ) + + yL[0:y.shape[0]] = y + return yL + else: + return y def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, verbose=True, x0=None, regularisation=None ) : @@ -852,6 +882,7 @@ cdef class Evaluation : nF = self.DICTIONARY['IC']['nF'] nE = self.DICTIONARY['EC']['nE'] nV = self.DICTIONARY['nV'] + nS = self.KERNELS['iso'].shape[1] norm_fib = np.ones( nF ) # x is the x of the original problem # self.x is the x preconditioned @@ -884,7 +915,7 @@ cdef class Evaluation : niiMAP_hdr['descrip'] = 'Created with COMMIT %s'%self.get_config('version') y_mea = np.reshape( self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float32), (nV,-1) ) - y_est = np.reshape( self.A.dot(self.x), (nV,-1) ).astype(np.float32) + y_est = np.reshape( self.A.dot(self.x)[0:int(nV*nS)], (nV,-1) ).astype(np.float32) print( '\t\t- RMSE... ', end='' ) sys.stdout.flush() diff --git a/commit/operator/operator.pyx b/commit/operator/operator.pyx index 6d83202a..abc4a784 100755 --- a/commit/operator/operator.pyx +++ b/commit/operator/operator.pyx @@ -26,6 +26,45 @@ cdef extern void COMMIT_At( unsigned char *_ICthreadsT, unsigned int *_ECthreadsT, unsigned int *_ISOthreadsT ) nogil +cdef extern void Tikhonov_L1( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L2( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L1z( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L2z( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L1t( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L2t( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L1zt( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil + +cdef extern void Tikhonov_L2zt( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_v_in, double *_v_out +) nogil cdef class LinearOperator : @@ -35,6 +74,8 @@ cdef class LinearOperator : """ cdef int nS, nF, nR, nE, nT, nV, nI, n, ndirs cdef public int adjoint, n1, n2 + cdef public float tikhonov_lambda + cdef public tikhonov_matrix cdef DICTIONARY cdef KERNELS @@ -61,20 +102,22 @@ cdef class LinearOperator : cdef unsigned int* ISOthreadsT - def __init__( self, DICTIONARY, KERNELS, THREADS ) : + def __init__( self, DICTIONARY, KERNELS, THREADS, tikhonov_lambda=0, tikhonov_matrix=None ) : """Set the pointers to the data structures used by the C code.""" self.DICTIONARY = DICTIONARY self.KERNELS = KERNELS self.THREADS = THREADS - self.nF = DICTIONARY['IC']['nF'] # number of FIBERS - self.nR = KERNELS['wmr'].shape[0] # number of FIBER RADII - self.nE = DICTIONARY['EC']['nE'] # number of EC segments - self.nT = KERNELS['wmh'].shape[0] # number of EC TORTUOSITY values - self.nV = DICTIONARY['nV'] # number of VOXELS - self.nI = KERNELS['iso'].shape[0] # number of ISO contributions - self.n = DICTIONARY['IC']['n'] # numbner of IC segments - self.ndirs = KERNELS['wmr'].shape[1] # number of directions + self.nF = DICTIONARY['IC']['nF'] # number of FIBERS + self.nR = KERNELS['wmr'].shape[0] # number of FIBER RADII + self.nE = DICTIONARY['EC']['nE'] # number of EC segments + self.nT = KERNELS['wmh'].shape[0] # number of EC TORTUOSITY values + self.nV = DICTIONARY['nV'] # number of VOXELS + self.nI = KERNELS['iso'].shape[0] # number of ISO contributions + self.n = DICTIONARY['IC']['n'] # numbner of IC segments + self.ndirs = KERNELS['wmr'].shape[1] # number of directions + self.tikhonov_lambda = tikhonov_lambda # equalizer parameter of the Tikhonov regularization term + self.tikhonov_matrix = tikhonov_matrix # derivative matrix of the Tikhonov regularization term if KERNELS['wmr'].size > 0 : self.nS = KERNELS['wmr'].shape[2] # number of SAMPLES @@ -85,7 +128,18 @@ cdef class LinearOperator : self.adjoint = 0 # direct of inverse product - self.n1 = self.nV*self.nS + # set shape of the operator according to tikhonov_matrix + if self.tikhonov_lambda > 0: + if self.tikhonov_matrix == 'L1': + self.n1 = self.nV*self.nS + (self.nR-1) + elif self.tikhonov_matrix == 'L2': + self.n1 = self.nV*self.nS + (self.nR-2) + elif self.tikhonov_matrix == 'L1z': + self.n1 = self.nV*self.nS + (self.nR+1) + elif self.tikhonov_matrix == 'L2z': + self.n1 = self.nV*self.nS + (self.nR) + else: + self.n1 = self.nV*self.nS self.n2 = self.nR*self.nF + self.nT*self.nE + self.nI*self.nV # get C pointers to arrays in DICTIONARY @@ -131,7 +185,7 @@ cdef class LinearOperator : @property def T( self ) : """Transpose of the explicit matrix.""" - C = LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS ) + C = LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, self.tikhonov_lambda, self.tikhonov_matrix ) C.adjoint = 1 - C.adjoint return C @@ -188,4 +242,58 @@ cdef class LinearOperator : self.ICthreadsT, self.ECthreadsT, self.ISOthreadsT ) + if self.tikhonov_lambda > 0: + if not self.adjoint: + # DIRECT PRODUCT lambda*L*x + if self.tikhonov_matrix == 'L1': + with nogil: + Tikhonov_L1( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L2': + with nogil: + Tikhonov_L2( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L1z': + with nogil: + Tikhonov_L1z( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L2z': + with nogil: + Tikhonov_L2z( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + else: + # INVERSE PRODUCT lambda*L'*y + if self.tikhonov_matrix == 'L1': + with nogil: + Tikhonov_L1t( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L2': + with nogil: + Tikhonov_L2t( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L1z': + with nogil: + Tikhonov_L1zt( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + elif self.tikhonov_matrix == 'L2z': + with nogil: + Tikhonov_L2zt( + self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda, + &v_in[0], &v_out[0] + ) + return v_out diff --git a/commit/operator/operator.pyxbld b/commit/operator/operator.pyxbld old mode 100644 new mode 100755 diff --git a/commit/operator/operator_noLUT.c b/commit/operator/operator_noLUT.c index 0046f237..84803cff 100644 --- a/commit/operator/operator_noLUT.c +++ b/commit/operator/operator_noLUT.c @@ -185,3 +185,5 @@ void COMMIT_At( pthread_join( threads[t], NULL ); return; } + +//TODO: Add tikhonov regularization when no LUT is required \ No newline at end of file diff --git a/commit/operator/operator_withLUT.c b/commit/operator/operator_withLUT.c index 9c959f57..effe2f2e 100644 --- a/commit/operator/operator_withLUT.c +++ b/commit/operator/operator_withLUT.c @@ -2245,3 +2245,139 @@ void COMMIT_At( pthread_join( threads[t], NULL ); return; } + +// =============================================================================================== +// Compute L*x MATRIX-VECTOR product (L is first order derivative with free boundary conditions) +// =============================================================================================== +void Tikhonov_L1( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int r = 0; r < _nIC-1; r++){ + for(int f = 0; f < _nF; f++){ + _vOUT[_nV*_nS + r] += _lambda*( -_vIN[r*_nF + f] + _vIN[(r+1)*_nF + f] ); + } + } +} + +// ========================================================================================================== +// Compute Lt*y TRANSPOSE-MATRIX-VECTOR product (L is first order derivative with free boundary conditions) +// ========================================================================================================== +void Tikhonov_L1t( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + _vOUT[f] += _lambda*( -_vIN[_nV*_nS] ); + + for(int r = 1; r < _nIC-1; r++) + _vOUT[_nF*r + f] += _lambda*( _vIN[_nV*_nS + r-1] - _vIN[_nV*_nS + r] ); + + _vOUT[_nF*(_nIC-1) + f] += _lambda*( _vIN[_nV*_nS + _nIC-2] ); + } +} + +// =============================================================================================== +// Compute L*x MATRIX-VECTOR product (L is second order derivative with free boundary conditions) +// =============================================================================================== +void Tikhonov_L2( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int r = 0; r < _nIC-2; r++){ + for(int f = 0; f < _nF; f++){ + _vOUT[_nV*_nS + r] += _lambda*( _vIN[r*_nF + f] -2*_vIN[(r+1)*_nF + f] + _vIN[(r+2)*_nF + f] ); + } + } +} + +// ========================================================================================================== +// Compute Lt*y TRANSPOSE-MATRIX-VECTOR product (L is second order derivative with free boundary conditions) +// ========================================================================================================== +void Tikhonov_L2t( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + _vOUT[f] += _lambda*( _vIN[_nV*_nS] ); + + _vOUT[_nF + f] += _lambda*( -2*_vIN[_nV*_nS] + _vIN[_nV*_nS + 1] ); + + for (int r = 2; r < _nIC-2; r++){ + _vOUT[r*_nF + f] += _lambda*( _vIN[_nV*_nS + (r-2)] -2*_vIN[_nV*_nS + (r-1)] + _vIN[_nV*_nS + r] ); + } + + _vOUT[(_nIC-2)*_nF + f] += _lambda*( _vIN[_nV*_nS + _nIC-4] -2*_vIN[_nV*_nS + _nIC-3] ); + + _vOUT[(_nIC-1)*_nF + f] += _lambda*( _vIN[_nV*_nS + (_nIC-3)] ); + } +} + +// =============================================================================================== +// Compute L*x MATRIX-VECTOR product (L is first order derivative with zero boundary conditions) +// =============================================================================================== +void Tikhonov_L1z( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + _vOUT[_nV*_nS] += _lambda*( _vIN[f] ); + + for(int r = 1; r < _nIC; r++){ + _vOUT[_nV*_nS + r] += _lambda*( -_vIN[(r-1)*_nF + f] + _vIN[r*_nF + f] ); + } + + _vOUT[_nV*_nS + _nIC] += _lambda*( -_vIN[(_nIC-1)*_nF + f] ); + } +} + +// ========================================================================================================== +// Compute Lt*y TRANSPOSE-MATRIX-VECTOR product (L is first order derivative with zero boundary conditions) +// ========================================================================================================== +void Tikhonov_L1zt( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + for(int r = 0; r < _nIC; r++){ + _vOUT[r*_nF + f] += _lambda*( _vIN[_nV*_nS + r] - _vIN[_nV*_nS + r + 1]); + } + } +} + +// =============================================================================================== +// Compute L*x MATRIX-VECTOR product (L is second order derivative with zero boundary conditions) +// =============================================================================================== +void Tikhonov_L2z( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + + _vOUT[_nV*_nS] += _lambda*( -2*_vIN[f] + x[_nF + f] ); + + for(int r = 1; r < _nIC-1; r++){ + _vOUT[_nV*_nS + r] += _lambda*( _vIN[(r-1)*_nF + f] -2*_vIN[r*_nF + f] + _vIN[(r+1)*_nF + f] ); + } + + _vOUT[_nV*_nS + _nIC - 1] += _lambda*( _vIN[(_nIC-2)*_nF + f] - 2*_vIN[(_nIC-1)*_nF + f] ); + } +} + +// ========================================================================================================== +// Compute Lt*y TRANSPOSE-MATRIX-VECTOR product (L is second order derivative with zero boundary conditions) +// ========================================================================================================== +void Tikhonov_L2zt( + int _nF, int _nIC, int _nV, int _nS, double _lambda, + double *_vIN, double *_vOUT) +{ + for(int f = 0; f < _nF; f++){ + _vOUT[f] += _lambda*( -2*_vIN[_nV*_nS] + _vIN[_nV*_nS + 1] ); + + for (int r = 0; r < _nIC-1; r++){ + _vOUT[r*_nF + f] += _lambda*( _vIN[_nV*_nS + (r-1)] - 2*_vIN[_nV*_nS + r] + _vIN[_nV*_nS + (r+1)] ); + } + + _vOUT[(_nIC-1)*_nF + f] += _lambda*( _vIN[_nV*_nS + (_nIC-2)] - 2*_vIN[_nV*_nS + (_nIC-1)] ); + } +} diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index d23e72fe..85f9f5c0 100755 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/commit/trk2dictionary/trk2dictionary.pyx @@ -426,4 +426,4 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam nii_hdr['descrip'] = 'Created with COMMIT %s'%get_distribution('dmri-commit').version nibabel.save( niiMASK, join(path_out,'dictionary_mask.nii.gz') ) - LOG( '\n [ %.1f seconds ]' % ( time.time() - tic ) ) \ No newline at end of file + LOG( '\n [ %.1f seconds ]' % ( time.time() - tic ) ) diff --git a/setup.py b/setup.py index 74391f2e..d8fabce4 100644 --- a/setup.py +++ b/setup.py @@ -44,7 +44,7 @@ def run(self): description = 'Convex Optimization Modeling for Microstructure Informed Tractography (COMMIT)' opts = dict(name='dmri-commit', - version='1.4.5', + version='1.5.0', description=description, long_description=description, author='Alessandro Daducci',