Skip to content

Commit

Permalink
Improve documentation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lmperry
lmperry committed Feb 7, 2017
1 parent 6b257be commit ea8c430
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 32 deletions.
39 changes: 37 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -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
```
30 changes: 13 additions & 17 deletions src/dtiError.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
%
Expand All @@ -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')
Expand All @@ -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);

Expand All @@ -84,25 +84,25 @@
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
bmask = niftiRead(wmProb);
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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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




31 changes: 18 additions & 13 deletions src/gear_dtiError.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...');

Expand Down Expand Up @@ -57,53 +62,53 @@ 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'])
ylabel([upper(eType{ii}), ' Predicted'])
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');
axis('image');axis('vis3d');
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);
Expand All @@ -122,4 +127,4 @@ function gear_dtiError(input_directory, output_directory, error_type, ncoords, w
exit(0)
end

return
return

0 comments on commit ea8c430

Please sign in to comment.