Message Boards Message Boards

Convolution Neural Network for MRI Video Data NDSB2 using NetChain?

Posted 8 years ago

Can NetChain and associated layers support 3 dimensional video data (MRI scan) as input or only really support 2D inputs (i.e. images)?

I am trying to apply the new Mathematica 11 deep neural network tools to process cardiac MRI slices from the second annual National Science Data Bowl.

The data is provided as a collection of "studies" (i.e. a unique patient heart) where each study has multiple scans each ~ 30x256x192 in size. That is, each patient heart is scanned multiple times from different angles/positions using MRI. Therefore each scan is a short "video clip" of a single heartbeat 30 frames long (i.e. about 1 second at 30 fps).

The ultimate goal is to predict the ratio of the minimum and maximum volume of blood being pumped in and out of the heart (i.e. (Vmax-Vmin)/Vmax). That is, what percentage of the blood into the heart is pumped out of the heart. Blood volume ratio too high or too low is indicative of a serious heart problem.

There are about 11000 total MRI scans (each 30x256x192 in size) with associated Vmin, Vmax labels.

I preprocessed all the DICOM files for the scans to normalize and resize the scans into 30x32x32.

Just trying to get something going, I used LeNet from the Mathematica docs with Vmin labels

net = NetChain[
  {
 ConvolutionLayer[20, {5, 5}],
 ElementwiseLayer[Ramp],
 PoolingLayer[{2, 2}, {2, 2}],
 ConvolutionLayer[50, {5, 5}],
 ElementwiseLayer[Ramp],
 PoolingLayer[{2, 2}, {2, 2}],
 FlattenLayer[],
 DotPlusLayer[500],
 ElementwiseLayer[Ramp],
 DotPlusLayer[1]       
},
  "Input" -> {30, 32, 32}, "Output" -> "Scalar"
]

However, I'm not sure that the dimensions are correct on the layers.

For example, when I just create a convolution layer

ConvolutionLayer[20,{5,5},"Input"->{30,32,32}]

The result says that the weights are a tensor of dimension {20,30,5,5} which seems right but then the biases are only of dimension {20} rather than {30,20} (?).

Then the ConvolutionLayer output is of dimension {20,28,28} but shouldn't the output be of dimension {30,20,28,28} ?

At any rate, the network trains fine (and fast with GPU as target device) but does not generalize at all (i.e. performs very poorly on hold out test scans). I'm wondering if the network actually has the appropriate configuration to really process the 30x32x32 input or if it's confused and only processing the first of the 30 images or something like that. Or maybe everything is just fine and that's just the way dimensions are being reported. I'm just guessing ...

Probably should rather be using softmax layer to predict "class" of the volume. That is, just bin the volume and predict the bin rather than have the network reproduce the numerical/decimal value for the volume.

Any help/advice/direction would be very much appreciated. Happy to share notebook and ResourceObject if anyone is interested.

Thanks

POSTED BY: Abel Brown
5 Replies
Posted 8 years ago

Thanks Daniel. This is awesome and very much appreciated. I think the performance of the predict functions might just roughly follow the distribution of the data (i.e. majority of the data has EF ~ 60% +- 5%). But, the method from your SYNASC paper is definitely really great. If interested, analysis of the top NDSB2 solutions (all used CNNs in workflow) here. Thanks again!

POSTED BY: Abel Brown

It turns out I can improve somewhat on the values shown in my response using Predict. For this I use a variant of a method presented at SYNASC 2016, " Linking Fourier and PCA Methods for Image Look-up".

http://synasc.ro/2016/list-of-papers/

The paper should appear in a few months. For now I think I should defer to that in terms of detailed explanations. But here is the rough idea.

(1) For each data item, take a discrete cosing transform. Keep only some low dimensional components. (2) Flatten the retained DCT matrix so each data item is now a vector. (3) Use the singular values decomposition on the matrix of flattened dcts that comes from the training set of data items. Keep only the largest singular values and corresponding singular vector matrices. (4) Use the result to create a lookup function. Apply it to test data elements to find some number of "closest" training vectors to each test element. Use the values of those nearby training elements to guess the test values.

The actual code is fairly simple and quite fast.

 nearestImages[dctvecs_, vals_, keep_] :=
 Module[
  {uu, ww, vv, udotw, norms},
  {uu, ww, vv} = SingularValueDecomposition[dctvecs, keep];
  udotw = uu.ww;
  norms = Map[Sqrt[#.#] &, udotw];
  udotw = udotw/norms;
  Clear[uu, ww];
  {Nearest[udotw -> vals(*,Method\[Rule]"KDTree"*)], vv}]

processInput[dctvecs_, vv_] :=
 Module[
  {tdotv, norms},
  tdotv = dctvecs.vv;
  norms = Map[Sqrt[#.#] &, tdotv];
  tdotv = tdotv/norms;
  tdotv]

guesses[nf_, tvecs_, n_] := Module[
  {nbrs, rng = Range[N@n, 1., -1]^2},
  nbrs = Map[nf[#, n] &, tvecs];
  Map[rng.# &, nbrs]/Total[rng]
  ]

Since this is fast, at least with the parameter settings I found to work well (based on smaller test runs), I use 10000 images for training.

n = 10000;
traindata = dat[[1 ;; n, 1]];
trainvalues = dat[[1 ;; n, 2]];
m = 1000;
testdata = dat[[n + 1 ;; n + m, 1]];
testvalues = dat[[n + 1 ;; n + m, 2]];
keep = 4;
dn = 4;
dct = 2;

AbsoluteTiming[
 traindctvecs = 
  Table[Flatten[
    Map[FourierDCT[#, dct] &, traindata[[j]]][[All, 1 ;; dn, 1 ;; dn]]],
   {j, Length[traindata]}];
 {nfunc, vv} =
  nearestImages[traindctvecs, trainvalues, keep];]
AbsoluteTiming[
 testdctvecs = 
  Table[Flatten[
    Map[FourierDCT[#, dct] &, testdata[[j]]][[All, 1 ;; dn, 1 ;; dn]]],
   {j, Length[testdata]}];
 testvecs = processInput[testdctvecs, vv];]
newvalsNI = guesses[nfunc, testvecs, 140];
relerrorsNI = 
  Abs[(newvalsNI - testvalues)/Sqrt[(newvalsNI^2 + testvalues)/2]];
{Length[Select[relerrorsNI, # > .2 &]],
  Length[Select[relerrorsNI, # > .4 &]],
  Length[Select[relerrorsNI, # < .1 &]]}/N[m]

(* Out[403]= {29.446372, Null}

Out[404]= {1.980268, Null}

Out[407]= {0.11, 0.031, 0.681} *)

So now we have 11% that are off by more than 20% of the correct values, 68% are within a tenth of the correct values, and only 3% have errors in excess of 40% of corresponding reference values. And it was several times faster than any of the Predict methods, when they trained on 2000 data values. Not bad for the end of the working day.

Since the speed is reasonable, one could of course do a number of randomized tests using RandomSample to separate the input into training vs test data. That could help to give some idea of the variance of this method, that is to say, how far it might stray from the values shown above.

POSTED BY: Daniel Lichtblau

I am not all that familiar with machine learning innards so I cannot help with specifics of NetChain. I did run some more basic code and will show the results here. First thing: there is a typograhic mismatch between your notebook and the reop directory layout. Very hard on those of us "of an age", in terms of trying to make sense of the error messages. Anyway, use "NDSB2" consistently and all will be better.

That said, the git clone instructions were excellent and very much appreciated. Especially by those of us "of an age", where figuring out how to properly install software quickly becomes hopeless.

I will pick up from where your code in ndsb2.nb has already set up dat. I use the first 2000 entries for training purposes and the next 1000 for testing. One could of course do many random runs since there are nearly 12000 entries. So here we proceed.

n = 2000;
traindata = dat[[1 ;; n, 1]];
trainvalues = dat[[1 ;; n, 2]];
m = 1000;
testdata = dat[[n + 1 ;; n + m, 1]];
testvalues = dat[[n + 1 ;; n + m, 2]];

I use several methods for the function Predict, and in each case use the setting PerformanceGoal -> "Quality". I show the timings both for setting up the predictor function, and for running it using each method setting. I also evaluate relative errors between predicted and actual values for the test entries. I give a list of the form {fraction worse than 20% off, fraction worse than 40% off, and fraction within 10%}. So we want values in the range 0 to 1 that we hope are of the form {small, very small, large}. I use the methods {"LinearRegression", "RandomForest", "NearestNeighbors", "NeuralNetwork"}. There is also a `"GaussianProcess" but it spewed errors when I tried it.

AbsoluteTiming[
 predictorLR = 
   Predict[traindata -> trainvalues, Method -> "LinearRegression", 
    PerformanceGoal -> "Quality"];]

(* Out[175]= {70.920143, Null} *)

AbsoluteTiming[newvalsLR = Map[predictorLR, testdata];]
relerrorsLR = 
  Abs[(newvalsLR - testvalues)/Sqrt[(newvalsLR^2 + testvalues)/2]];
{Length[Select[relerrorsLR, # > .2 &]],
  Length[Select[relerrorsLR, # > .4 &]],
  Length[Select[relerrorsLR, # < .1 &]]}/N[m]

(* Out[191]= {11.857996, Null}

Out[193]= {0.229, 0.072, 0.501} *)

AbsoluteTiming[
 predictorRF = 
   Predict[traindata -> trainvalues, Method -> "RandomForest", 
    PerformanceGoal -> "Quality"];]

(* Out[179]= {163.405545, Null} *)

AbsoluteTiming[newvalsRF = Map[predictorRF, testdata];]
relerrorsRF = 
  Abs[(newvalsRF - testvalues)/Sqrt[(newvalsRF^2 + testvalues)/2]];
{Length[Select[relerrorsRF, # > .2 &]],
  Length[Select[relerrorsRF, # > .4 &]],
  Length[Select[relerrorsRF, # < .1 &]]}/N[m]

(* Out[194]= {13.969906, Null}

Out[196]= {0.195, 0.077, 0.542} *)

AbsoluteTiming[
 predictorNN = 
   Predict[traindata -> trainvalues, Method -> "NearestNeighbors", 
    PerformanceGoal -> "Quality"];]

(* Out[183]= {68.806587, Null} *)

AbsoluteTiming[newvalsNN = Map[predictorNN, testdata];]
relerrorsNN = 
  Abs[(newvalsNN - testvalues)/Sqrt[(newvalsNN^2 + testvalues)/2]];
{Length[Select[relerrorsNN, # > .2 &]],
  Length[Select[relerrorsNN, # > .4 &]],
  Length[Select[relerrorsNN, # < .1 &]]}/N[m]

(* Out[197]= {12.105393, Null}

Out[199]= {0.371, 0.169, 0.366} *)

AbsoluteTiming[
 predictorNN2 = 
   Predict[traindata -> trainvalues, Method -> "NeuralNetwork", 
    PerformanceGoal -> "Quality"];]

(* Out[187]= {455.508036, Null} *)

AbsoluteTiming[newvalsNN2 = Map[predictorNN2, testdata];]
relerrorsNN2 = 
  Abs[(newvalsNN2 - testvalues)/Sqrt[(newvalsNN2^2 + testvalues)/2]];
{Length[Select[relerrorsNN2, # > .2 &]],
  Length[Select[relerrorsNN2, # > .4 &]],
  Length[Select[relerrorsNN2, # < .1 &]]}/N[m]

(* Out[200]= {22.463613, Null}

Out[202]= {0.231, 0.073, 0.505} *)

Here is a summary of the results.

(1) "NearestNeighbors" gives a poor result relative to the rest. We discard it from further consideration.

(2) "NeuralNetwork" method is slow.

(3) All three of the better performers have around 7-8% that are worse than 40% off the mark, 20-23% that are worse than 20% off, and 50-54% that come within 10% of correct values.

This strikes me as being a reasonable outcome though again, I do not have enough familiarity with either the methods or this data to say whether one might expect to do better.

POSTED BY: Daniel Lichtblau
Posted 8 years ago

Yes, repo is here on git hub (~800 MB). Do a git clone

git clone https://github.com/softwarespartan/NDSB2.git

and then run the notebook ndsb2.nb. You'll need to set the "root" variable in the notebook to the local directory with the repository.

POSTED BY: Abel Brown

Do you have a repository with the 30x32x32 scans? Not necessarily all 11000 of them, but maybe 1000 or so? (If all are available somewhere that would be fine too.)

POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract