Merge pull request #117 from OSIPI/matlab_GU

All algorithms wrapped and brain test data

Author: oliverchampion

.github/workflows/unit_test.yml

@@ -35,4 +35,5 @@ jobs:
         run: |
           pip install pytest pytest-cov
           python -m pytest --doctest-modules --junitxml=junit/test-results.xml --cov=. --cov-report=xml --cov-report=html -r w
+          python -m pytest --doctest-modules --junitxml=junit/test-results-brain.xml --cov=. --cov-report=xml --cov-report=html --dataFile tests/IVIMmodels/unit_tests/generic_brain.json -r w -k test_ivim_fit_saved
 

.gitignore

@@ -5,6 +5,8 @@
 __pycache__/
 *.nii.gz
 *.nii
+*.bval
+*.bvec
 *.npy
 *.dcm
 *.npy
@@ -37,3 +39,4 @@ coverage.xml
 phantoms/MR_XCAT_qMRI/*.json
 phantoms/MR_XCAT_qMRI/*.txt
 tests/IVIMmodels/unit_tests/models
+models

README.md

@@ -22,6 +22,9 @@ If you would like to contribute with code, please follow the instructions below:
 If you would like to use code from the repository and/or are new to Github or IVIM, please see the jupyter notebook below:
 *   [Introduction to TF2.4_IVIM-MRI_CodeCollection github and IVIM Analysis using Python](doc/Introduction_to_TF24_IVIM-MRI_CodeCollection_github_and_IVIM_Analysis_using_Python.ipynb)
 
+If you would like to use MATLAB-based algorithms or apply the testing framework to one of these algorithms, please have a look at these instructions:
+* [MATLAB-related instructions](doc/matlab_instructions.md)
+
 ## Repository Organization
 
 The repository is organized in four main folders along with configuration files for automated testing. 

conftest.py

@@ -264,8 +264,7 @@ def algorithmlist(algorithms):
     for algorithm in algorithms["algorithms"]:
         algorithm_dict = algorithms.get(algorithm, {})
         requires_matlab = algorithm_dict.get("requires_matlab", False)
-        if not algorithm_dict.get('deep_learning', False):
-            yield algorithm, requires_matlab
+        yield algorithm, requires_matlab, algorithm_dict.get('deep_learning', False)
 
 
 def bound_input(datafile, algorithms):

doc/matlab_instructions.md

@@ -0,0 +1,7 @@
+## MATLAB-related instructions
+
+A large portion of IVIM-related code has historically been written in the MATLAB programming language. For this reason, we have included the possibility to call and test MATLAB-based code from the Python-based code developed here. This requires a few additional steps to get going
+
+- In addition to the Python packages listed requirements.txt, one also need to install matlab.engine, for example with pip: `python -m pip install matlabengine`. This requires a corresponding version of MATLAB to be installed.
+- To be able to call the MATLAB code in this repository, the src/original folder must be added to the MATLAB path
+- Testing of MATLAB-based code only runs offline (not on Github). This is done by running `python -m pytest --withmatlab`

phantoms/brain/data/diffusive_snr2000.bval

@@ -0,0 +1 @@
+0 10 20 30 40 80 110 140 170 200 300 400 500 600 700 800 900

requirements.txt

@@ -17,4 +17,5 @@ sphinx
 sphinx_rtd_theme
 pytest-json-report
 statsmodels
-ivimnet
+ivimnet
+nlopt

src/original/OJ_GU/IVIM_bayes.m

@@ -181,7 +181,7 @@
 for i = 1:ceil(M/voxelsPerRun)
     % partition voxels
     usedVoxels = (i-1)*voxelsPerRun + 1:min(i*voxelsPerRun,M);
-    
+
     fpart = f(usedVoxels);
     Dpart = D(usedVoxels);
     Dstarpart = Dstar(usedVoxels);
@@ -263,11 +263,11 @@
 
         % Display iteration every 500th iteration
         if ~mod(j,500) && j > burns
-            disp(['Iterations: ' num2str(j-burns)]);
+            %disp(['Iterations: ' num2str(j-burns)]);
         elseif ~mod(j,100) && j < burns
-            disp(['Burn in-steps: ' num2str(j)]);
+            %disp(['Burn in-steps: ' num2str(j)]);
         elseif j == burns
-            disp(['Burn in complete: ' num2str(j)]);
+            %disp(['Burn in complete: ' num2str(j)]);
         end
     end
 
@@ -286,28 +286,28 @@
         out.S0.std(usedVoxels) = sqrt(theta2sum(:,4)/n-(thetasum(:,4)/n).^2);
     else
         %mean
-        out.f.mean(usedVoxels) = mean(squeeze(theta(:,1,burns + 1:n+burns)),2);
-        out.D.mean(usedVoxels) = mean(squeeze(theta(:,2,burns + 1:n+burns)),2);
-        out.Dstar.mean(usedVoxels) = mean(squeeze(theta(:,3,burns + 1:n+burns)),2);
-        out.S0.mean(usedVoxels) = mean(squeeze(theta(:,4,burns + 1:n+burns)),2);
+        out.f.mean(usedVoxels) = mean(reshape(theta(:,1,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.D.mean(usedVoxels) = mean(reshape(theta(:,2,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.Dstar.mean(usedVoxels) = mean(reshape(theta(:,3,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.S0.mean(usedVoxels) = mean(reshape(theta(:,4,burns + 1:n+burns),[size(theta,1),n]),2);
 
         %median
-        out.f.median(usedVoxels) = median(squeeze(theta(:,1,burns + 1:n+burns)),2);
-        out.D.median(usedVoxels) = median(squeeze(theta(:,2,burns + 1:n+burns)),2);
-        out.Dstar.median(usedVoxels) = median(squeeze(theta(:,3,burns + 1:n+burns)),2);
-        out.S0.median(usedVoxels) = median(squeeze(theta(:,4,burns + 1:n+burns)),2);
+        out.f.median(usedVoxels) = median(reshape(theta(:,1,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.D.median(usedVoxels) = median(reshape(theta(:,2,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.Dstar.median(usedVoxels) = median(reshape(theta(:,3,burns + 1:n+burns),[size(theta,1),n]),2);
+        out.S0.median(usedVoxels) = median(reshape(theta(:,4,burns + 1:n+burns),[size(theta,1),n]),2);
 
         % mode
-        out.f.mode(usedVoxels) = halfSampleMode(squeeze(theta(:,1,burns + 1:n+burns)));
-        out.D.mode(usedVoxels) = halfSampleMode(squeeze(theta(:,2,burns + 1:n+burns)));
-        out.Dstar.mode(usedVoxels) = halfSampleMode(squeeze(theta(:,3,burns + 1:n+burns)));
-        out.S0.mode(usedVoxels) = halfSampleMode(squeeze(theta(:,4,burns + 1:n+burns)));
+        out.f.mode(usedVoxels) = halfSampleMode(reshape(theta(:,1,burns + 1:n+burns),[size(theta,1),n]));
+        out.D.mode(usedVoxels) = halfSampleMode(reshape(theta(:,2,burns + 1:n+burns),[size(theta,1),n]));
+        out.Dstar.mode(usedVoxels) = halfSampleMode(reshape(theta(:,3,burns + 1:n+burns),[size(theta,1),n]));
+        out.S0.mode(usedVoxels) = halfSampleMode(reshape(theta(:,4,burns + 1:n+burns),[size(theta,1),n]));
 
         % standard deviation
-        out.f.std(usedVoxels) = std(squeeze(theta(:,1,burns + 1:n+burns)),1,2);
-        out.D.std(usedVoxels) = std(squeeze(theta(:,2,burns + 1:n+burns)),1,2);
-        out.Dstar.std(usedVoxels) = std(squeeze(theta(:,3,burns + 1:n+burns)),1,2);
-        out.S0.std(usedVoxels) = std(squeeze(theta(:,4,burns + 1:n+burns)),1,2);
+        out.f.std(usedVoxels) = std(reshape(theta(:,1,burns + 1:n+burns),[size(theta,1),n]),1,2);
+        out.D.std(usedVoxels) = std(reshape(theta(:,2,burns + 1:n+burns),[size(theta,1),n]),1,2);
+        out.Dstar.std(usedVoxels) = std(reshape(theta(:,3,burns + 1:n+burns),[size(theta,1),n]),1,2);
+        out.S0.std(usedVoxels) = std(reshape(theta(:,4,burns + 1:n+burns),[size(theta,1),n]),1,2);
     end
 end
 
@@ -381,7 +381,52 @@
 y(b < 3.75) = log(1.0 + a1.*(3.5156229 + a1.*(3.0899424 + a1.*(1.2067492 + a1.*(0.2659732 + a1.*(0.0360768 + a1.*0.0045813))))));
 y(b >= 3.75) = b(b >= 3.75) + log((0.39894228+a2.*(0.01328592+a2.*(0.00225319+a2.*(-0.00157565+a2.*(0.00916281+a2.*(-0.02057706+a2.*(0.02635537+a2.*(-0.01647633+a2.*0.00392377))))))))./sqrt(b(b>=3.75)));
 
-
+function hsm = halfSampleMode(X)
+% calculates the half sample mode for each row in X
+if ~ismatrix(X)
+    error('X must be a matrix or vector');
+end
+X = sort(X,2);
+n = size(X,2);
+hsm = HSM_rec(n,X);
+function hsm = HSM_rec(n,X)
+% special cases
+if size(X,2) == 1
+    hsm = X;
+    return;
+elseif size(X,2) == 2
+    hsm = sum(X,2)/2;
+    return;
+elseif size(X,2) == 3   
+    hsm = zeros(size(X,1),1);
+    low = (X(:,2) - X(:,1)) < (X(:,3) - X(:,2)); % use lower pair
+    eq = (X(:,2) - X(:,1)) == (X(:,3) - X(:,2)); % use mid value
+    
+    if any(low)
+        hsm(low) = sum(X(low,1:2),2)/2;
+    end
+    if any(eq)
+        hsm(eq) = X(eq,2);
+    end
+    if any(~(low|eq))
+        hsm(~(low|eq)) = sum(X(~(low|eq),2:3),2)/2; % otherwise
+    end
+    return
+end
+% general case (n > 3)
+wmin = X(:,end) - X(:,1);
+N = ceil(n/2);
+j = ones(size(X,1),1);
+for i = 1:(n-N+1)
+   w = X(:,i+N-1) - X(:,i);
+   m = w < wmin;
+   wmin(m) = w(m);
+   j(m) = i;
+end
+rowsub = repmat((1:size(X,1))',1,N);
+colsub = repmat(j,1,N) + repmat(0:N-1,size(X,1),1);
+Xsub = reshape(X(sub2ind(size(X),rowsub,colsub)),size(X,1),N);
+hsm = HSM_rec(N,Xsub);
 
 
 

src/original/OJ_GU/IVIM_seg.m

@@ -140,12 +140,13 @@
     A = sum(Yred.*exp(-bred*D))./sum(exp(-2*bred*D));
 
     % Assures that A is within limits
-    D(A < lim(1,2)) = lim(1,1);
-    D(A > lim(2,2)) = lim(2,1);
-    maskD(A < lim(1,2) | A > lim(2,2)) = false;
-    A(A < lim(1,2)) = lim(1,2);
-    A(A > lim(2,2)) = lim(2,2);
-    
+    %D(A < lim(1,2)) = lim(1,1);
+    %D(A > lim(2,2)) = lim(2,1);
+    %maskD(A < lim(1,2) | A > lim(2,2)) = false;
+    A(A < 0) = 0;
+    %A(A < lim(1,2)) = lim(1,2);
+    %A(A > lim(2,2)) = lim(2,2);
+
     % Calculate f
     if estf && ~step2
         if disp_prog
@@ -173,31 +174,36 @@
 
         % calculate possible range of f to remove voxels where f cannot
         % result in a fit within limits
-        f1 = zeros(1,n);
-        f1(mask_back) = fcalc(lim(1,4)*ones(1,sum(mask_back)),mask_back);
-        f2 = zeros(1,n);
-        f2(mask_back) = fcalc(lim(2,4)*ones(1,sum(mask_back)),mask_back);
+        %f1 = zeros(1,n);
+        %f1(mask_back) = fcalc(lim(1,4)*ones(1,sum(mask_back)),mask_back);
+        %f2 = zeros(1,n);
+        %f2(mask_back) = fcalc(lim(2,4)*ones(1,sum(mask_back)),mask_back);
 
-        maskf = maskD & (((f1 > lim(1,3)) & (f1 < lim(2,3))) | ((f2 > lim(1,3)) & (f2 < lim(2,3))));
-        if disp_prog
-            fprintf('Discarding %d voxels due to f out of bounds\n',sum(maskD) - sum(maskf));
-        end
+        %maskf = maskD & (((f1 > lim(1,3)) & (f1 < lim(2,3))) | ((f2 > lim(1,3)) & (f2 < lim(2,3))));
+        maskf = mask_back;
+        %if disp_prog
+        %    fprintf('Discarding %d voxels due to f out of bounds\n',sum(maskD) - sum(maskf));
+        %end
 
         % Prepares output
         Dstar = zeros(1,n); 
-        Dstar(f1 < lim(1,3)) = lim(1,4);
-        Dstar(f2 > lim(2,3)) = lim(2,4);
+        %Dstar(f1 < lim(1,3)) = lim(1,4);
+        %Dstar(f2 > lim(2,3)) = lim(2,4);
         maskDstar = maskf;
 
         % Optimization
-        [Dstar(maskf),maskDstar(maskf)] = optimizeD(Y(:,maskf)-repmat(A(maskf),length(b),1).*exp(-b*D(maskf)),b,optlim,lim(:,4),disp_prog);
-    
+        if sum(maskf) > 0
+            [Dstar(maskf),maskDstar(maskf)] = optimizeD(Y(:,maskf)-repmat(A(maskf),length(b),1).*exp(-b*D(maskf)),b,optlim,lim(:,4),disp_prog);
+        end
+        
         % Calculates f
         f = lim(1,1)*ones(1,n);
-        f(f1 < lim(1,3)) = lim(1,3);
-        f(f2 > lim(2,3)) = lim(2,3); 
-        f(maskDstar) = fcalc(Dstar(maskDstar),maskDstar);
-    
+        %f(f1 < lim(1,3)) = lim(1,3);
+        %f(f2 > lim(2,3)) = lim(2,3); 
+        if sum(maskDstar) > 0
+            f(maskDstar) = fcalc(Dstar(maskDstar),maskDstar);
+        end
+
         % Checks for f out of bounds
         maskf = maskf & (f > lim(1,3)) & (f < lim(2,3));
         Dstar(f < lim(1,3)) = lim(1,4);

src/standardized/IAR_LU_biexp.py

@@ -79,7 +79,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
 
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
-            gtab = gradient_table(bvalues, bvec, b0_threshold=0)
+            gtab = gradient_table(bvalues, bvecs=bvec, b0_threshold=0)
 
             self.IAR_algorithm = IvimModelBiExp(gtab, bounds=self.bounds, initial_guess=self.initial_guess)
 
@@ -111,7 +111,7 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
 
             bvec = np.zeros((bvalues.size, 3))
             bvec[:,2] = 1
-            gtab = gradient_table(bvalues, bvec, b0_threshold=0)
+            gtab = gradient_table(bvalues, bvecs=bvec, b0_threshold=0)
 
             self.IAR_algorithm = IvimModelBiExp(gtab, bounds=self.bounds, initial_guess=self.initial_guess)
         b0_index = np.where(bvalues == 0)[0][0]