2D CNN to classify 3D Greyscale MRI data, possible data labeling issue

Question

I am trying to run a binary classification of 3D black and white MRI data. I'm using 2D convolutions due to the lack of channels inherent in B&W data. I've added a dimension to get the dimensionality to line up, and essentially the depth of this data is acting as the batch dimension. I am using a subsample of data, 20 files each 189 by 233 by 197. Just as a quick background.

I have a csv file with a bunch of info including the label data for each file which I've tried to extract, as per the code below.

import numpy as np
import glob
import os
import tensorflow as tf
import pandas as pd
import glob

import SimpleITK as sitk

from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from keras.preprocessing.image import ImageDataGenerator


from keras.utils import plot_model
from keras.utils import to_categorical
from keras.utils import np_utils

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout

from google.colab import drive
drive.mount('/content/gdrive')

datapath = ('/content/gdrive/My Drive/DirectoryTest/All Data/')
patients = os.listdir(datapath)
labels_df = pd.read_csv('/content/Data_Index.csv', index_col = 0 )

labelset = []

for i in patients:
  label = labels_df.loc[i, 'Group']
  if label is 'AD':
    np.char.replace(label, ['AD'], [0])
  if label is 'CN':
    np.char.replace(label, ['CN'], [1])
  labelset.append(label)

label_encoder = LabelEncoder()
labelset = label_encoder.fit_transform(labelset)

labelset = np_utils.to_categorical(labelset, num_classes= 2)

FullDataSet = []

for i in patients:
  a = sitk.ReadImage(datapath + i)
  b = sitk.GetArrayFromImage(a)
  c = np.reshape(b, (189,233,197, 1))
  FullDataSet.append(c)

training_data, testing_data, training_labels, testing_labels = train_test_split(FullDataSet, labelset, train_size=0.70,test_size=0.30)

dataset_train = tf.data.Dataset.from_tensor_slices((training_data, training_labels))
dataset_test = tf.data.Dataset.from_tensor_slices((testing_data, testing_labels))

CNN_model = tf.keras.Sequential(
  [
      #tf.keras.layers.Input(shape=(189, 233, 197, 1), batch_size=2),
      #tf.keras.layers.Reshape((197, 233, 189, 1)),   
                              
      tf.keras.layers.Conv2D(kernel_size=(7, 7), data_format='channels_last', filters=64, activation='relu',
                             padding='same', strides=( 3, 3), input_shape=( 233, 197, 1)),
      #tf.keras.layers.BatchNormalization(center=True, scale=False),
      tf.keras.layers.MaxPool2D(pool_size=(3, 3), padding='same'),
      tf.keras.layers.Dropout(0.20),
      
      tf.keras.layers.Conv2D(kernel_size=( 7, 7), filters=128, activation='relu', padding='same', strides=( 3, 3)),
      #tf.keras.layers.BatchNormalization(center=True, scale=False),
      tf.keras.layers.MaxPool2D(pool_size=(3, 3), padding='same'),
      tf.keras.layers.Dropout(0.20),      

      tf.keras.layers.Conv2D(kernel_size=( 7, 7), filters=256, activation='relu', padding='same', strides=( 3, 3)),
      #tf.keras.layers.BatchNormalization(center=True, scale=False),
      tf.keras.layers.MaxPool2D(pool_size=(3, 3), padding = 'same'),
      tf.keras.layers.Dropout(0.20), 

      # last activation could be either sigmoid or softmax, need to look into this more. Sig for binary output, Soft for multi output 
      tf.keras.layers.Flatten(),
      tf.keras.layers.Dense(256, activation='relu'),   
      tf.keras.layers.Dense(64, activation='relu'),
      tf.keras.layers.Dropout(0.20),
      tf.keras.layers.Dense(2, activation='softmax')

  ])
# Compile the model
CNN_model.compile(optimizer=tf.keras.optimizers.Adam(lr=0.00001), loss='binary_crossentropy', metrics=['accuracy'])

# print model layers
CNN_model.summary()

CNN_history = CNN_model.fit(dataset_train, epochs=10, validation_data=dataset_test)

When I go to fit the model I get the following error:

Epoch 1/10
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-a8b210ec2e72> in <module>()
      1 #running of the model
      2 #CNN_history = CNN_model.fit(dataset_train, epochs=100, validation_data =dataset_test, validation_steps=1)
----> 3 CNN_history = CNN_model.fit(dataset_train, epochs=10, validation_data=dataset_test)
      4 
      5 

10 frames
/usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/func_graph.py in wrapper(*args, **kwargs)
    971           except Exception as e:  # pylint:disable=broad-except
    972             if hasattr(e, "ag_error_metadata"):
--> 973               raise e.ag_error_metadata.to_exception(e)
    974             else:
    975               raise

ValueError: in user code:

    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/engine/training.py:806 train_function  *
        return step_function(self, iterator)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/engine/training.py:796 step_function  **
        outputs = model.distribute_strategy.run(run_step, args=(data,))
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/distribute/distribute_lib.py:1211 run
        return self._extended.call_for_each_replica(fn, args=args, kwargs=kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/distribute/distribute_lib.py:2585 call_for_each_replica
        return self._call_for_each_replica(fn, args, kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/distribute/distribute_lib.py:2945 _call_for_each_replica
        return fn(*args, **kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/engine/training.py:789 run_step  **
        outputs = model.train_step(data)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/engine/training.py:749 train_step
        y, y_pred, sample_weight, regularization_losses=self.losses)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/engine/compile_utils.py:204 __call__
        loss_value = loss_obj(y_t, y_p, sample_weight=sw)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/losses.py:149 __call__
        losses = ag_call(y_true, y_pred)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/losses.py:253 call  **
        return ag_fn(y_true, y_pred, **self._fn_kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/util/dispatch.py:201 wrapper
        return target(*args, **kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/losses.py:1605 binary_crossentropy
        K.binary_crossentropy(y_true, y_pred, from_logits=from_logits), axis=-1)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/util/dispatch.py:201 wrapper
        return target(*args, **kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/keras/backend.py:4829 binary_crossentropy
        bce = target * math_ops.log(output + epsilon())
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/math_ops.py:1141 binary_op_wrapper
        raise e
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/math_ops.py:1125 binary_op_wrapper
        return func(x, y, name=name)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/math_ops.py:1457 _mul_dispatch
        return multiply(x, y, name=name)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/util/dispatch.py:201 wrapper
        return target(*args, **kwargs)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/math_ops.py:509 multiply
        return gen_math_ops.mul(x, y, name)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/gen_math_ops.py:6176 mul
        "Mul", x=x, y=y, name=name)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/op_def_library.py:744 _apply_op_helper
        attrs=attr_protos, op_def=op_def)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/func_graph.py:593 _create_op_internal
        compute_device)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/ops.py:3485 _create_op_internal
        op_def=op_def)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/ops.py:1975 __init__
        control_input_ops, op_def)
    /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/ops.py:1815 _create_c_op
        raise ValueError(str(e))

    ValueError: Dimensions must be equal, but are 2 and 189 for '{{node binary_crossentropy/mul}} = Mul[T=DT_FLOAT](ExpandDims, binary_crossentropy/Log)' with input shapes: [2,1], [189,2].

I know the 2 in [189,2] is associated with the final softmax layer but I don't know what to do with that information, or where to go from here. Any help would be appreciated, thanks!

Answer 1

Here are some comments regarding your code, hopefully helpful.

Use Conv3D and MaxPool3D

If you're dealing with 3D images, then you should almost certainly use Conv3D instead of Conv2D, and MaxPool3D instead of MaxPool2D. Here's an example (using random data), which I just tested, and it seems to work fine:

import numpy as np
import tensorflow as tf
from tensorflow import keras

train_size = 20
val_size = 5

X_train = np.random.random([train_size, 189, 233, 197]).astype(np.float32)
X_valid = np.random.random([val_size, 189, 233, 197]).astype(np.float32)
y_train = np.random.randint(2, size=train_size).astype(np.float32)
y_valid = np.random.randint(2, size=val_size).astype(np.float32)

CNN_model = keras.Sequential([
      keras.layers.Reshape([189, 233, 197, 1], input_shape=[189, 233, 197]),
      keras.layers.Conv3D(kernel_size=(7, 7, 7), filters=32, activation='relu',
                          padding='same', strides=(3, 3, 3)),
      #keras.layers.BatchNormalization(),
      keras.layers.MaxPool3D(pool_size=(3, 3, 3), padding='same'),
      keras.layers.Dropout(0.20),
      
      keras.layers.Conv3D(kernel_size=(5, 5, 5), filters=64, activation='relu',
                          padding='same', strides=(3, 3, 3)),
      #keras.layers.BatchNormalization(),
      keras.layers.MaxPool3D(pool_size=(2, 2, 2), padding='same'),
      keras.layers.Dropout(0.20),

      keras.layers.Conv3D(kernel_size=(3, 3, 3), filters=128, activation='relu',
                          padding='same', strides=(1, 1, 1)),
      #keras.layers.BatchNormalization(),
      keras.layers.MaxPool3D(pool_size=(2, 2, 2), padding='same'),
      keras.layers.Dropout(0.20),

      keras.layers.Flatten(),
      keras.layers.Dense(256, activation='relu'),   
      keras.layers.Dense(64, activation='relu'),
      keras.layers.Dropout(0.20),
      keras.layers.Dense(1, activation='sigmoid')
  ])

# Compile the model
CNN_model.compile(optimizer=keras.optimizers.Adam(lr=0.00001),
                  loss='binary_crossentropy',
                  metrics=['accuracy'])

# print model layers
CNN_model.summary()

CNN_history = CNN_model.fit(X_train, y_train, epochs=10,
                            validation_data=[X_valid, y_valid])

Don't reshape to permute dimensions

Regarding these two commented out lines:

      #tf.keras.layers.Input(shape=(189, 233, 197, 1), batch_size=2),
      #tf.keras.layers.Reshape((197, 233, 189, 1)),   

Reshaping a 189x233x197x1 image to 197x233x189x1 will not work as you expect. It will completely shuffle the pixels around, making the task much harder. This is similar to reshaping a 2x3 image into a 3x2 image:

>>> img = np.array([[1,2,3],[4,5,6]])
>>> np.reshape(img, [3, 2])
array([[1, 2],
       [3, 4],
       [5, 6]])

Notice that this is not the same as rotating the image: the pixels are completely mixed up.

What you want is to use tf.keras.layers.Permute() instead like this:

CNN_model = tf.keras.Sequential([
      tf.keras.layers.Permute((3, 2, 1, 4), input_shape=(189, 233, 197, 1)),
      ...
])

Because these commented out lines were wrong, I suspect that the following line may also be wrong:

c = np.reshape(b, (189,233,197, 1))

I don't know the shape of b, so please make absolutely sure that it is compatible with this np.reshape() operation. For example, if its shape is [189, 233, 197], that's fine. But if it is [197, 233, 189], for example, then you will need to permute dimensions before you reshape:

b_permuted = np.transpose(b, [2, 1, 0]) # permute dims
c = np.reshape(b_permuted, [189, 233, 197, 1]) # then add the channels dim

The np.transpose() function is similar to using Permute(), except dimensions are 0-indexed instead of 1-indexed.

It may be even more complex. For example, if the 3D images are stored as large 2D images containing smaller 2D slices side by side, then the shape of b may be something like [189*197, 233]. In this case, you would need to do something like this:

b_reshaped = np.reshape(b, [189, 197, 233, 1])
c = np.transpose(b_reshaped, [0, 2, 1, 3])

I hope these examples are clear enough.

Use tf.keras, not keras

There are several implementations of the Keras API. One is the keras package, which is the "multi-backend" Keras (which is installed using pip install keras). Another is tf.keras and it comes with TensorFlow. Your program seems to use both. You should absolutely avoid that, it will cause weird issues.

from keras.utils import plot_model # this is multibackend Keras
...
CNN_model = tf.keras.Sequential(...) # this is tf.keras

I highly recommend you uninstall multibackend keras, to avoid this type of error: pip uninstall keras. Then fix the imports by prefixing with tensorflow., for example:

from tensorflow.keras.models import Sequential
from tensorflow.keras.utils import to_categorical # note: not from np_utils
...

Don't use to_categorical() for binary classification

For binary classification, the labels should simply be a 1D array containing 0s and 1s, such as np.array([1., 0., 0., 1., 1.]). The code can be quite simplified:

labelset = []

for i in patients:
  label = labels_df.loc[i, 'Group']
  if label == 'AD':  # use `==` instead of `is` to compare strings
    labelset.append(0.)
  elif label == 'CN':
    labelset.append(1.)
  else:
      raise "Oops, unknown label" # I recommend testing possible failure cases

labelset = np.array(labelset)

Importantly, for binary classification, you should use a single neuron in the output layer, and also use the "sigmoid" activation function (not "softmax", which is used for multiclass classification):

CNN_model = tf.keras.Sequential([
      ...
      tf.keras.layers.Dense(1, activation='sigmoid')
])

Minor comment

• You don't need to specify both the train_size and the test_size when calling train_test_split().

Good luck!