annotate linear_fascile_evaluation.py @ 1:84a2e30b5404 draft

Uploaded
author greg
date Tue, 28 Nov 2017 13:30:57 -0500
parents cbfa8c336751
children 0ddfcb3b5ce6
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
cbfa8c336751 Uploaded
greg
parents:
diff changeset
1 #!/usr/bin/env python
cbfa8c336751 Uploaded
greg
parents:
diff changeset
2 import argparse
cbfa8c336751 Uploaded
greg
parents:
diff changeset
3 import numpy as np
cbfa8c336751 Uploaded
greg
parents:
diff changeset
4 import os.path as op
cbfa8c336751 Uploaded
greg
parents:
diff changeset
5 import nibabel as nib
cbfa8c336751 Uploaded
greg
parents:
diff changeset
6 import dipy.core.optimize as opt
cbfa8c336751 Uploaded
greg
parents:
diff changeset
7 import dipy.tracking.life as life
cbfa8c336751 Uploaded
greg
parents:
diff changeset
8 import matplotlib.pyplot as plt
cbfa8c336751 Uploaded
greg
parents:
diff changeset
9 import matplotlib
cbfa8c336751 Uploaded
greg
parents:
diff changeset
10
cbfa8c336751 Uploaded
greg
parents:
diff changeset
11 from dipy.viz.colormap import line_colors
cbfa8c336751 Uploaded
greg
parents:
diff changeset
12 from dipy.viz import fvtk
cbfa8c336751 Uploaded
greg
parents:
diff changeset
13 from mpl_toolkits.axes_grid1 import AxesGrid
1
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
14 from dipy.data import read_stanford_labels, fetch_stanford_t1, read_stanford_t1
0
cbfa8c336751 Uploaded
greg
parents:
diff changeset
15
cbfa8c336751 Uploaded
greg
parents:
diff changeset
16 parser = argparse.ArgumentParser()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
17 parser.add_argument('--candidates', dest='candidates', help='Candidates selection')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
18 parser.add_argument('--output_life_candidates', dest='output_life_candidates', help='Output life candidates')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
19 parser.add_argument('--output_life_optimized', dest='output_life_optimized', help='Output life optimized streamlines')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
20 parser.add_argument('--output_beta_histogram', dest='output_beta_histogram', help='Output beta histogram')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
21 parser.add_argument('--output_error_histograms', dest='output_error_histograms', help='Output error histograms')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
22 parser.add_argument('--output_spatial_errors', dest='output_spatial_errors', help='Output spatial errors')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
23
cbfa8c336751 Uploaded
greg
parents:
diff changeset
24 args = parser.parse_args()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
25
1
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
26 # We'll need to know where the corpus callosum is from these variables.
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
27 hardi_img, gtab, labels_img = read_stanford_labels()
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
28 labels = labels_img.get_data()
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
29 cc_slice = labels == 2
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
30 fetch_stanford_t1()
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
31 t1 = read_stanford_t1()
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
32 t1_data = t1.get_data()
84a2e30b5404 Uploaded
greg
parents: 0
diff changeset
33 data = hardi_img.get_data()
0
cbfa8c336751 Uploaded
greg
parents:
diff changeset
34
cbfa8c336751 Uploaded
greg
parents:
diff changeset
35 # Read the candidates from file in voxel space:
cbfa8c336751 Uploaded
greg
parents:
diff changeset
36 candidate_sl = [s[0] for s in nib.trackvis.read(args.candidates, points_space='voxel')[0]]
cbfa8c336751 Uploaded
greg
parents:
diff changeset
37 # Visualize the initial candidate group of streamlines
cbfa8c336751 Uploaded
greg
parents:
diff changeset
38 # in 3D, relative to the anatomical structure of this brain.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
39 candidate_streamlines_actor = fvtk.streamtube(candidate_sl, line_colors(candidate_sl))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
40 cc_ROI_actor = fvtk.contour(cc_slice, levels=[1], colors=[(1., 1., 0.)], opacities=[1.])
cbfa8c336751 Uploaded
greg
parents:
diff changeset
41 vol_actor = fvtk.slicer(t1_data)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
42 vol_actor.display(40, None, None)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
43 vol_actor2 = vol_actor.copy()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
44 vol_actor2.display(None, None, 35)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
45 # Add display objects to canvas.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
46 ren = fvtk.ren()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
47 fvtk.add(ren, candidate_streamlines_actor)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
48 fvtk.add(ren, cc_ROI_actor)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
49 fvtk.add(ren, vol_actor)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
50 fvtk.add(ren, vol_actor2)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
51 fvtk.record(ren, n_frames=1, out_path=args.output_life_candidates, size=(800, 800))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
52 # Initialize a LiFE model.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
53 fiber_model = life.FiberModel(gtab)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
54 # Fit the model, producing a FiberFit class instance,
cbfa8c336751 Uploaded
greg
parents:
diff changeset
55 # that stores the data, as well as the results of the
cbfa8c336751 Uploaded
greg
parents:
diff changeset
56 # fitting procedure.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
57 fiber_fit = fiber_model.fit(data, candidate_sl, affine=np.eye(4))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
58 fig, ax = plt.subplots(1)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
59 ax.hist(fiber_fit.beta, bins=100, histtype='step')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
60 ax.set_xlabel('Fiber weights')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
61 ax.set_ylabel('# fibers')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
62 fig.savefig(args.output_beta_histogram)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
63 # Filter out these redundant streamlines and
cbfa8c336751 Uploaded
greg
parents:
diff changeset
64 # generate an optimized group of streamlines.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
65 optimized_sl = list(np.array(candidate_sl)[np.where(fiber_fit.beta>0)[0]])
cbfa8c336751 Uploaded
greg
parents:
diff changeset
66 ren = fvtk.ren()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
67 fvtk.add(ren, fvtk.streamtube(optimized_sl, line_colors(optimized_sl)))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
68 fvtk.add(ren, cc_ROI_actor)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
69 fvtk.add(ren, vol_actor)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
70 fvtk.record(ren, n_frames=1, out_path=args.output_life_optimized, size=(800, 800))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
71 model_predict = fiber_fit.predict()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
72 # Focus on the error in prediction of the diffusion-weighted
cbfa8c336751 Uploaded
greg
parents:
diff changeset
73 # data, and calculate the root of the mean squared error.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
74 model_error = model_predict - fiber_fit.data
cbfa8c336751 Uploaded
greg
parents:
diff changeset
75 model_rmse = np.sqrt(np.mean(model_error[:, 10:] ** 2, -1))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
76 # Calculate another error term by assuming that the weight for each streamline
cbfa8c336751 Uploaded
greg
parents:
diff changeset
77 # is equal to zero. This produces the naive prediction of the mean of the
cbfa8c336751 Uploaded
greg
parents:
diff changeset
78 # signal in each voxel.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
79 beta_baseline = np.zeros(fiber_fit.beta.shape[0])
cbfa8c336751 Uploaded
greg
parents:
diff changeset
80 pred_weighted = np.reshape(opt.spdot(fiber_fit.life_matrix, beta_baseline), (fiber_fit.vox_coords.shape[0], np.sum(~gtab.b0s_mask)))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
81 mean_pred = np.empty((fiber_fit.vox_coords.shape[0], gtab.bvals.shape[0]))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
82 S0 = fiber_fit.b0_signal
cbfa8c336751 Uploaded
greg
parents:
diff changeset
83 # Since the fitting is done in the demeaned S/S0 domain,
cbfa8c336751 Uploaded
greg
parents:
diff changeset
84 # add back the mean and then multiply by S0 in every voxel:
cbfa8c336751 Uploaded
greg
parents:
diff changeset
85 mean_pred[..., gtab.b0s_mask] = S0[:, None]
cbfa8c336751 Uploaded
greg
parents:
diff changeset
86 mean_pred[..., ~gtab.b0s_mask] = (pred_weighted + fiber_fit.mean_signal[:, None]) * S0[:, None]
cbfa8c336751 Uploaded
greg
parents:
diff changeset
87 mean_error = mean_pred - fiber_fit.data
cbfa8c336751 Uploaded
greg
parents:
diff changeset
88 mean_rmse = np.sqrt(np.mean(mean_error ** 2, -1))
cbfa8c336751 Uploaded
greg
parents:
diff changeset
89 # Compare the overall distribution of errors
cbfa8c336751 Uploaded
greg
parents:
diff changeset
90 # between these two alternative models of the ROI.
cbfa8c336751 Uploaded
greg
parents:
diff changeset
91 fig, ax = plt.subplots(1)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
92 ax.hist(mean_rmse - model_rmse, bins=100, histtype='step')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
93 ax.text(0.2, 0.9,'Median RMSE, mean model: %.2f' % np.median(mean_rmse), horizontalalignment='left', verticalalignment='center', transform=ax.transAxes)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
94 ax.text(0.2, 0.8,'Median RMSE, LiFE: %.2f' % np.median(model_rmse), horizontalalignment='left', verticalalignment='center', transform=ax.transAxes)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
95 ax.set_xlabel('RMS Error')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
96 ax.set_ylabel('# voxels')
cbfa8c336751 Uploaded
greg
parents:
diff changeset
97 fig.savefig(args.output_error_histograms)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
98 # Show the spatial distribution of the two error terms,
cbfa8c336751 Uploaded
greg
parents:
diff changeset
99 # and of the improvement with the model fit:
cbfa8c336751 Uploaded
greg
parents:
diff changeset
100 vol_model = np.ones(data.shape[:3]) * np.nan
cbfa8c336751 Uploaded
greg
parents:
diff changeset
101 vol_model[fiber_fit.vox_coords[:, 0], fiber_fit.vox_coords[:, 1], fiber_fit.vox_coords[:, 2]] = model_rmse
cbfa8c336751 Uploaded
greg
parents:
diff changeset
102 vol_mean = np.ones(data.shape[:3]) * np.nan
cbfa8c336751 Uploaded
greg
parents:
diff changeset
103 vol_mean[fiber_fit.vox_coords[:, 0], fiber_fit.vox_coords[:, 1], fiber_fit.vox_coords[:, 2]] = mean_rmse
cbfa8c336751 Uploaded
greg
parents:
diff changeset
104 vol_improve = np.ones(data.shape[:3]) * np.nan
cbfa8c336751 Uploaded
greg
parents:
diff changeset
105 vol_improve[fiber_fit.vox_coords[:, 0], fiber_fit.vox_coords[:, 1], fiber_fit.vox_coords[:, 2]] = mean_rmse - model_rmse
cbfa8c336751 Uploaded
greg
parents:
diff changeset
106 sl_idx = 49
cbfa8c336751 Uploaded
greg
parents:
diff changeset
107 fig = plt.figure()
cbfa8c336751 Uploaded
greg
parents:
diff changeset
108 fig.subplots_adjust(left=0.05, right=0.95)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
109 ax = AxesGrid(fig, 111, nrows_ncols = (1, 3), label_mode = "1", share_all = True, cbar_location="top", cbar_mode="each", cbar_size="10%", cbar_pad="5%")
cbfa8c336751 Uploaded
greg
parents:
diff changeset
110 ax[0].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
111 im = ax[0].matshow(np.rot90(vol_model[sl_idx, :, :]), cmap=matplotlib.cm.hot)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
112 ax.cbar_axes[0].colorbar(im)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
113 ax[1].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
114 im = ax[1].matshow(np.rot90(vol_mean[sl_idx, :, :]), cmap=matplotlib.cm.hot)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
115 ax.cbar_axes[1].colorbar(im)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
116 ax[2].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
117 im = ax[2].matshow(np.rot90(vol_improve[sl_idx, :, :]), cmap=matplotlib.cm.RdBu)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
118 ax.cbar_axes[2].colorbar(im)
cbfa8c336751 Uploaded
greg
parents:
diff changeset
119 for lax in ax:
cbfa8c336751 Uploaded
greg
parents:
diff changeset
120 lax.set_xticks([])
cbfa8c336751 Uploaded
greg
parents:
diff changeset
121 lax.set_yticks([])
cbfa8c336751 Uploaded
greg
parents:
diff changeset
122 fig.savefig(args.output_spatial_errors)