From ea8c43007cdce4ad468427c155ee3c0f77d7b977 Mon Sep 17 00:00:00 2001 From: Michael Perry Date: Tue, 7 Feb 2017 12:14:46 -0800 Subject: [PATCH] Improve documentation. --- README.md | 39 +++++++++++++++++++++++++++++++++++++-- src/dtiError.m | 30 +++++++++++++----------------- src/gear_dtiError.m | 31 ++++++++++++++++++------------- 3 files changed, 68 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index 5481f47..843340c 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,37 @@ -# scitran/dtiError -Calculate the RMSE between a tensor fit from dtiInit and the diffusion weighted imaging data +# scitran/dti-error +Calculate RMSE between the measured signal and the ADC (or dSIG) based on tensor model fit provided by [dtiInit](https://github.com/scitran-apps/dtiinit). + +This gear calculates the histogram of differences between DTI based predictions (ADC or dSig) with the actual ADC or dSig data. Larger deviations suggest noisier data. This is one of a series of methods we are developing to assess the reliability of diffusion weighted image data. + +### Build the Gear +```#bash +git clone https://github.com/scitran-apps/dti-error +cd dti-error +./build.sh +``` + +### Example Usage +First download a [dtiInit](https://github.com/scitran-apps/dtiinit) output archive and save to disk (e.g., dtiInit_27-Jan-2017_18-51-24.zip used below). Then run the commands below: + +```bash +# Directory and file handling +DTIINIT_OUTPUT=dtiInit_27-Jan-2017_18-51-24.zip # Use the name of your dtiInit output +INPUT_DIR=`pwd`/input/dtiInit_Archive +OUTPUT_DIR=`pwd`/output +mkdir -p $INPUT_DIR +mkdir -p $OUTPUT_DIR + +# Run the Gear +docker run --rm -ti -v $INPUT_DIR:/flywheel/v0/input -v $OUTPUT_DIR:/flywheel/v0/output scitran/dti-error:v0.1.0 +``` + +### Compile the Matlab Executable +Clone required code and prepare for compilation: +```bash +git clone https://github.com/scitran-apps/dti-error +git clone https://github.com/vistalab/vistasoft +``` +In Matlab (e.g., r2015b): +```Matlab +mcc -m dti-error/src/gear_dtiError.m -I vistasoft -I dti-error/src +``` diff --git a/src/dtiError.m b/src/dtiError.m index 2148ef8..8e47fa9 100644 --- a/src/dtiError.m +++ b/src/dtiError.m @@ -32,12 +32,12 @@ % % [X Y Z] = meshgrid(40, 40, 40); coords = [X(:) Y(:) Z(:)]; % % err = dtiError(baseName,'coords',coords,'eType','adc'); -% mrvNewGraphWin; +% mrvNewGraphWin; % hist(err,50); xlabel('\Delta ADC'); ylabel('Count') % fprintf('DWI image quality %.2f (ADC-DTI method, higher better)\n',1/std(err)); % % err = dtiError(baseName,'coords',coords,'eType','dsig'); -% mrvNewGraphWin; +% mrvNewGraphWin; % hist(err,50); xlabel('\Delta DSIG'); ylabel('Count') % fprintf('DWI image quality %.2f (DSIG-DTI eval method, higher better)\n',1/std(err)); % @@ -46,7 +46,7 @@ % mrvNewGraphWin; % hist(err,50); xlabel('\Delta ADC'); ylabel('Count') % fprintf('DWI image quality %.2f (ADC-DTI method, higher better)\n',1/std(err)); -% +% % err = dtiError(baseName,'wmProb',wmProb,'eType','dsig','ncoords',250); % mrvNewGraphWin; % hist(err,50); xlabel('\Delta DSIG'); ylabel('Count') @@ -69,7 +69,7 @@ wmProb = p.Results.wmProb; if ~isempty(wmProb) && ~exist(wmProb,'file') - error('White matter probability file %s not found',wmProb); + error('White matter probability file %s not found',wmProb); end % ni = niftiRead(wmProb); niftiView(ni); @@ -84,7 +84,7 @@ if isempty(coords) if ~isempty(wmProb) % If there is a brain mask passed in use ncoords within the brain - % mask. + % mask. ncoords = p.Results.ncoords; % We read and randomly define selected values here @@ -92,17 +92,17 @@ bmask.data(bmask.data <= 180) = 0; bmask.data(bmask.data > 180) = 1; % niftiView(bmask); - + [i,j,k] = ind2sub(size(bmask.data),find(bmask.data == 1)); imax = length(i); lst = randi(imax,ncoords,1); allCoords = [i,j,k]; coords = allCoords(lst,:); - + else % Not sure what to do here, yet. We might find the b=0 data, select % high values, and choose nCoords within that - disp('No coords and no wmProb. Returning'); + disp('No coords and no wmProb. Returning'); return; % Get the b=0 data and find the nCoords from the top 50% of the top % b values @@ -127,13 +127,13 @@ % Instead, we separately calculate the predicted ADC values bvecs = dwiGet(dwi,'diffusion bvecs'); adcPredicted = dtiADC(Q,bvecs); - + % A visualization, if you like % dwiPlot(dwi,'adc',adc(:,5),squeeze(Q(:,:,5))); % mrvNewGraphWin; % plot(adc(:),adcPredicted(:),'o') % identityLine(gca); axis equal - + % Sometimes data are drawn from just outside of the white matter % region and the ADC values are not valid. This happens from time % to time, and we aren't worrying about it. We just delete those @@ -151,7 +151,7 @@ % It is possible to return the predicted ADC from dwiQ, as well, % by a small adjustment to that function. Q = dwiQ(dwi,coords); - + % Sometimes there are bad coordinates out there. We protect you. nBadADCcoords = sum(isnan(Q(1,1,:))); if nBadADCcoords @@ -168,18 +168,14 @@ S0 = dwiGet(dwi,'b0valsimage',coords); dsigPredicted = dwiComputeSignal(S0, bvecs, bvals, Q); dsigPredicted = dsigPredicted'; - + % We could use percentage error by dividing by dsig(). err = dsig(:) - dsigPredicted(:); measured = dsig(:); predicted = dsigPredicted(:); - + otherwise error('Unknown error type %s\n',eType); end end - - - - diff --git a/src/gear_dtiError.m b/src/gear_dtiError.m index 0e1893e..129fcc1 100644 --- a/src/gear_dtiError.m +++ b/src/gear_dtiError.m @@ -15,8 +15,13 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w % gear_dtiError(input_directory, output_directory, error_type, ncoords, wm_prob) % % EXAMPLE COMPILATION: -% mcc -m gear_dtiError.m -I ~/Code/vistasoft +% git clone https://github.com/scitran-apps/dti-error +% git clone https://github.com/vistalab/vistasoft % +% In Matlab (e.g., r2015b): +% mcc -m dti-error/src/gear_dtiError.m -I vistasoft -I dti-error/src +% + disp('Alive...'); @@ -57,16 +62,16 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w else [err, dwi, coords, measured, predicted] = dtiError(diffusion_nifti, 'eType', eType{ii}, 'ncoords', ncoords); end - + % PLOTS - + % Histogram mrvNewGraphWin; hist(err,50); xlabel(['\Delta ', upper(eType{ii})]); ylabel('Voxel Count') fprintf('%s: DWI image quality (1/std(err)) = %.2f \n', upper(eType{ii}), 1/std(err)); title(sprintf('%s: DWI image quality (1/std(err)) = %.2f \n', upper(eType{ii}), 1/std(err))); saveas(gcf, fullfile(output_directory, [eType{ii}, '_', num2str(ncoords), '_err.png'])); - + % Scatter with id/line plot(measured, predicted,'o') xlabel([upper(eType{ii}), ' Measured']) @@ -74,18 +79,18 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w title([upper(eType{ii}), ' Measured/Predicted']); identityLine(gca); axis equal saveas(gcf, fullfile(output_directory, [eType{ii}, '_', num2str(ncoords), '_measured_predicted.png'])); - - + + switch lower(eType{ii}) case 'adc' % Plot surface and points adc = dwiGet(dwi,'adc data image',coords); Q = dwiQ(dwi,coords); dwiPlot(dwi,'adc',adc(:,5),squeeze(Q(:,:,5))); - + % Write out gif of surface vs measured gif_file = fullfile(output_directory, [ eType{ii}, '_', num2str(ncoords), '_surface.gif' ] ); - + % Turn of the axes and loop through the views lightH = camlight('right'); axis('off'); @@ -93,17 +98,17 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w for ff = 1:180 % rotate the camera camorbit(2, 0); - + % Rotate the light with the camera lightH = camlight(lightH,'right'); - + % Caputure the frame mov(ff) = getframe(gcf); - + % Convert the movie frame to an image im = frame2im(mov(ff)); [imind, cm] = rgb2ind(im, 256); - + % Write out the image to gif if ff == 1; imwrite(imind, cm, gif_file, 'gif', 'Loopcount',inf, 'delaytime',0); @@ -122,4 +127,4 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w exit(0) end -return \ No newline at end of file +return