corrected [x,y,z] order in train and target_generation. Makes more sense for...

corrected [x,y,z] order in train and target_generation. Makes more sense for user now. Still in 'wrong' order in segment, but has no incidence on user. Also finished comments in deepfind.py

README.md

@@ -4,13 +4,15 @@ The code in this repository is described in [this pre-print](https://hal.inria.f
 
 __Disclaimer:__ this is a preliminary version of the code, which is subject to improvements. For ease of use, we also plan to create a graphical interface in forthcoming months.
 
+__News__: (29/01/20) A first version of the GUI is now available in folder pyqt/. [More information...](###Using the GUI) 
+
 ## Contents
 - [System requirements](##System requirements)
 - [Installation guide](##Installation guide)
 - [Instructions for use](##Instructions for use)
 
 ## System requirements
-__Deep Finder__ has been implemented using __Python 2.7__ and is based on the __Keras__ package. It has been tested on Linux (Debian 8.6), and should also work on Mac OSX as well as Windows.
+__Deep Finder__ has been implemented using __Python 3__ and is based on the __Keras__ package. It has been tested on Linux (Debian 8.6), and should also work on Mac OSX as well as Windows.
 
 The algorithm needs an Nvidia GPU to run at reasonable speed (in particular for training). The present code has been tested on Tesla K80 and M40 GPUs. For running on other GPUs, some parameter values (e.g. patch and batch sizes) may need to be changed to adapt to available memory.
 
@@ -21,33 +23,47 @@ The algorithm needs an Nvidia GPU to run at reasonable speed (in particular for
 ### Package dependencies
 Users should install following packages in order to run Deep Finder. The package versions for which our software has been tested are displayed in brackets:
 ```
-tensorflow-gpu (1.4.0)
-keras          (2.1.6)
-numpy          (1.14.3)
-h5py           (2.7.1)
-lxml           (4.3.2)
-scikit-learn   (0.19.1)     
-scikit-image   (0.14.2)  
-matplotlib     (2.2.3)
+tensorflow-gpu (1.14.0)
+keras          (2.3.1)
+numpy          (1.16.4)
+h5py           (2.9.0)
+lxml           (4.3.4)
+scikit-learn   (0.21.2)     
+scikit-image   (0.15.0)  
+matplotlib     (3.1.0)
+mrcfile        (1.1.2)
+PyQt5          (5.13.2)
 ```
 
 ## Installation guide
 First install the packages:
 ```
-pip install numpy tensorflow-gpu keras sklearn h5py lxml scikit-learn scikit-image matplotlib
+pip install numpy tensorflow-gpu keras sklearn h5py lxml scikit-learn scikit-image matplotlib mrcfile PyQt5
 ```
 For more details about installing Keras, please see [Keras installation instructions](https://keras.io/#installation).
 
 Once the dependencies are installed, the user should be able to run Deep Finder.
 
 ## Instructions for use
-Instructions for using Deep Finder are contained in folder examples/. The scripts contain comments on how the toolbox should be used. To run a script, first launch ipython:
+### Using the scripts
+Instructions for using Deep Finder are contained in folder examples/. The scripts contain comments on how the toolbox should be used. To run a script, first place yourself in its folder. For example, to run the target generation script:
 ```
-ipython 
+cd examples/training/
+python step1_generate_target.py
 ```
-Then launch the script:
+
+### Using the GUI
+The GUI (Graphical User Interface) is available in folder pyqt/, and should be more intuitive for those who are not used to work with script. For now, 5 GUIs are available (target generation, training, segmentation, clustering) and allow the same functionalities as the scripts in example/. To run a GUI, first place yourself in its folder. For example, to run the target generation GUI:
 ```
-%run script_file.py
+cd pyqt/generate_target/
+python gui_target.py
 ```
 
-__Note:__ working examples are contained in examples/analyze/, where Deep Finder processes the test tomogram from the [SHREC'19 challenge](http://www2.projects.science.uu.nl/shrec/cryo-et/). The script in examples/training/ will fail because the training data is not included in this Gitlab. In addition, the evaluation script (examples/analyze/step3_launch_evaluation.py) is the one used in SHREC'19, which needs python3 and additional packages. The performance of Deep Finder has been evaluated by an independent group, and the result of this evaluation has been published in Gubins & al., "SHREC'19 track: Classification in cryo-electron tomograms".
+![Training GUI](./screenshots/gui_train.png)
+
+In near future we plan to provide a data visualization tool.
+
+__Notes:__ 
+- working examples are contained in examples/analyze/, where Deep Finder processes the test tomogram from the [SHREC'19 challenge](http://www2.projects.science.uu.nl/shrec/cryo-et/). 
+- The script in examples/training/ will fail because the training data is not included in this Gitlab. 
+- The evaluation script (examples/analyze/step3_launch_evaluation.py) is the one used in SHREC'19, which needs additional packages (pathlib and pycm, can be installed with pip). The performance of Deep Finder has been evaluated by an independent group, and the result of this evaluation has been published in Gubins & al., "SHREC'19 track: Classification in cryo-electron tomograms".

core_utils.py

@@ -11,10 +11,13 @@ import matplotlib.pyplot as plt
 # INPUTS:
 #   path_data  : list of strings '/path/to/tomogram.ext'
 #   path_target: list of strings '/path/to/target.ext'
+#                The idx of above lists correspond to each other so that (path_data[idx], path_target[idx]) corresponds
+#                to a (tomog, target) pair
 #   dset_name  : can be usefull if files are stored as .h5
 # OUTPUTS:
 #   data_list  : list of 3D numpy arrays (tomograms)
 #   target_list: list of 3D numpy arrays (annotated tomograms)
+#                In the same way as for the inputs, (data_list[idx],target_list[idx]) corresponds to a (tomo,target) pair
 def load_dataset(path_data, path_target, dset_name='dataset'):
     data_list   = []
     target_list = []
@@ -23,7 +26,14 @@ def load_dataset(path_data, path_target, dset_name='dataset'):
         target_list.append(utils.read_array(path_target[idx], dset_name))
     return data_list, target_list
 
-
+# This function applies bootstrap (i.e. re-sampling) in case of unbalanced classes.
+# Given an objlist containing objects from various classes, this function outputs an equal amount of objects for each
+# class, each objects being uniformely sampled inside its class set.
+# INPUTS:
+#   objlist: list of dictionaries
+#   Nbs    : number of objects to sample from each class
+# OUTPUT:
+#   bs_idx : list of indexes corresponding to the bootstraped objects
 def get_bootstrap_idx(objlist,Nbs):
     # Get a vector containing the object class labels (from objlist):
     Nobj = len(objlist)
@@ -37,10 +47,19 @@ def get_bootstrap_idx(objlist,Nbs):
     # ->from label_list, sample Nbs objects from each class
     bs_idx = []
     for l in lblTAB:
-        bs_idx.append( np.random.choice(np.squeeze(np.asarray(np.nonzero(label_list==l))), Nbs) )
+        bs_idx.append( np.random.choice(np.squeeze(np.asarray(np.nonzero(label_list==l))), Nbs) ) # TODO: can be simplified
     bs_idx = np.concatenate(bs_idx)
     return bs_idx
-    
+
+# Takes position specified in 'obj', applies random shift to it, and then checks if the patch around this position is
+# out of the tomogram boundaries. If so, the position is shifted to that patch is inside the tomo boundaries.
+# INPUTS:
+#   tomodim: tuple (dimX,dimY,dimZ) containing size of tomogram
+#   p_in   : int lenght of patch in voxels
+#   obj    : dictionary obtained when calling objlist[idx]
+#   Lrnd   : int random shift in voxels applied to position
+# OUTPUTS:
+#   x,y,z  : int,int,int coordinates for sampling patch safely
 def get_patch_position(tomodim, p_in, obj, Lrnd):
     # sample at coordinates specified in obj=objlist[idx]
     x = int( obj['x'] )
@@ -56,9 +75,9 @@ def get_patch_position(tomodim, p_in, obj, Lrnd):
     if (x<p_in) : x = p_in
     if (y<p_in) : y = p_in
     if (z<p_in) : z = p_in
-    if (x>tomodim[0]-p_in): x = tomodim[0]-p_in
+    if (x>tomodim[2]-p_in): x = tomodim[2]-p_in
     if (y>tomodim[1]-p_in): y = tomodim[1]-p_in
-    if (z>tomodim[2]-p_in): z = tomodim[2]-p_in
+    if (z>tomodim[0]-p_in): z = tomodim[0]-p_in
 
     #else: # sample random position in tomogram
     #    x = np.int32( np.random.choice(range(p_in,tomodim[0]-p_in)) )
@@ -66,7 +85,11 @@ def get_patch_position(tomodim, p_in, obj, Lrnd):
     #    z = np.int32( np.random.choice(range(p_in,tomodim[0]-p_in)) )
     
     return x,y,z
-    
+
+# Saves training history as .h5 file.
+# INPUTS:
+#   history: dictionary object containing lists. These lists contain scores and metrics wrt epochs.
+#   filename: string '/path/to/net_train_history.h5'
 def save_history(history, filename):
     h5file = h5py.File(filename, 'w')
 
@@ -90,7 +113,11 @@ def save_history(history, filename):
 
     h5file.close()
     return
-    
+
+# Plots the training history as several graphs and saves them in an image file.
+# INPUTS:
+#   history: dictionary object containing lists. These lists contain scores and metrics wrt epochs.
+#   filename: string '/path/to/net_train_history_plot.png'
 def plot_history(history, filename):
     Ncl = len(history['val_f1'][0])
     legend_names = []
@@ -136,7 +163,7 @@ def plot_history(history, filename):
     fig.savefig(filename)
 
 
-# Following observer classes are needed to send prints to GUI if needed:
+# Following observer classes are needed to send prints to GUI:
 class observer_print:
     def display(message):
         print(message)

deepfind.py

@@ -15,7 +15,13 @@ import utils
 import core_utils
 import utils_objl as ol
 
-# TODO: add path_output as class argument
+# Note:
+# All imports relative to keras&tensorflow are done inside the constructor of classes that need it (i.e. Train and
+# Segment). This is not best practice, but importing tensorflow-related utilities takes time, and I only want to import
+# them when needed.
+# -> does not work cause 'to_categorical' is used in train method. This method does not see the import in class
+# constructor
+
 class DeepFinder:
     def __init__(self):
         self.obs_list = [core_utils.observer_print]
@@ -46,6 +52,7 @@ class TargetBuilder(DeepFinder):
     # INPUTS
     #   objl: list of dictionaries. Needs to contain [phi,psi,the] Euler angles for orienting the shapes.
     #   target_array: 3D numpy array that initializes the training target. Allows to pass an array already containing annotated structures like membranes.
+    #                 index order of array should be [z,y,x]
     #   ref_list: list of binary 3D arrays (expected to be cubic). These reference arrays contain the shape of macromolecules ('1' for 'is object' and '0' for 'is not object')
     #             The references order in list should correspond to the class label
     #             For ex: 1st element of list -> reference of class 1
@@ -66,7 +73,7 @@ class TargetBuilder(DeepFinder):
             the = objl[p]['the']
 
             ref = ref_list[lbl - 1]
-            centeroffset = np.int(np.floor(ref.shape[0] / 2))
+            centeroffset = np.int(np.floor(ref.shape[0] / 2)) # here we expect ref to be cubic
 
             # Rotate ref:
             if phi!=None and psi!=None and the!=None:
@@ -75,16 +82,17 @@ class TargetBuilder(DeepFinder):
 
             # Get the coordinates of object voxels in target_array
             obj_voxels = np.nonzero(ref == 1)
-            x_vox = obj_voxels[0] + x - centeroffset +1
+            x_vox = obj_voxels[2] + x - centeroffset +1
             y_vox = obj_voxels[1] + y - centeroffset +1
-            z_vox = obj_voxels[2] + z - centeroffset +1
+            z_vox = obj_voxels[0] + z - centeroffset +1
 
             for idx in range(x_vox.size):
                 xx = x_vox[idx]
                 yy = y_vox[idx]
                 zz = z_vox[idx]
-                if xx >= 0 and xx < dim[0] and yy >= 0 and yy < dim[1] and zz >= 0 and zz < dim[2]:  # if in tomo bounds
-                    target_array[xx, yy, zz] = lbl
+                if xx >= 0 and xx < dim[2] and yy >= 0 and yy < dim[1] and zz >= 0 and zz < dim[0]:  # if in tomo bounds
+                    target_array[zz, yy, xx] = lbl
+
         return np.int8(target_array)
 
     # Generates segmentation targets from object list. Here macromolecules are annotated with spheres.
@@ -93,6 +101,7 @@ class TargetBuilder(DeepFinder):
     # INPUTS
     #   objl: list of dictionaries. Needs to contain [phi,psi,the] Euler angles for orienting the shapes.
     #   target_array: 3D numpy array that initializes the training target. Allows to pass an array already containing annotated structures like membranes.
+    #                 index order of array should be [z,y,x]
     #   radius_list: list of sphere radii (in voxels).
     #             The radii order in list should correspond to the class label
     #             For ex: 1st element of list -> sphere radius for class 1
@@ -131,20 +140,32 @@ class Train(DeepFinder):
         self.flag_batch_bootstrap = 0
         self.Lrnd = 13  # random shifts applied when sampling data- and target-patches (in voxels)
 
-    # This function launches the training procedure. For each epoch, an image is plotted, displaying the progression with different metrics: loss, accuracy, f1-score, recall, precision. Every 10 epochs, the current network weights are saved.
+        # Build network:
+        self.net.compile(optimizer=self.optimizer, loss=losses.tversky_loss, metrics=['accuracy'])
+
+    # This function launches the training procedure. For each epoch, an image is plotted, displaying the progression
+    # with different metrics: loss, accuracy, f1-score, recall, precision. Every 10 epochs, the current network weights
+    # are saved.
     # INPUTS:
     #   path_data     : a list containing the paths to data files (i.e. tomograms)
     #   path_target   : a list containing the paths to target files (i.e. annotated volumes)
-    #   objlist_train : an xml structure containing information about annotated objects: origin tomogram (should correspond to the index of 'path_data' argument), coordinates, class. During training, these coordinates are used for guiding the patch sampling procedure.
-    #   objlist_valid : same as 'objlist_train', but objects contained in this xml structure are not used for training, but for validation. It allows to monitor the training and check for over/under-fitting. Ideally, the validation objects should originate from different tomograms than training objects.
-    # The network is trained on small 3D patches (i.e. sub-volumes), sampled from the larger tomograms (due to memory limitation). The patch sampling is not realized randomly, but is guided by the macromolecule coordinates contained in so-called object lists (objlist).
+    #   objlist_train : list of dictionaries containing information about annotated objects (e.g. class, position)
+    #                   In particular, the tomo_idx should correspond to the index of 'path_data' and 'path_target'.
+    #                   See utils_objl.py for more info about object lists.
+    #                   During training, these coordinates are used for guiding the patch sampling procedure.
+    #   objlist_valid : same as 'objlist_train', but objects contained in this list are not used for training,
+    #                   but for validation. It allows to monitor the training and check for over/under-fitting. Ideally,
+    #                   the validation objects should originate from different tomograms than training objects.
+    # The network is trained on small 3D patches (i.e. sub-volumes), sampled from the larger tomograms (due to memory
+    # limitation). The patch sampling is not realized randomly, but is guided by the macromolecule coordinates contained
+    # in so-called object lists (objlist).
     # Concerning the loading of the dataset, two options are possible:
     #    flag_direct_read=0: the whole dataset is loaded into memory
-    #    flag_direct_read=1: only the patches are loaded into memory, each time a training batch is generated. This is usefull when the dataset is too large to load into memory. However, the transfer speed between the data server and the GPU host should be high enough, else the procedure becomes very slow.
+    #    flag_direct_read=1: only the patches are loaded into memory, each time a training batch is generated. This is
+    #                        usefull when the dataset is too large to load into memory. However, the transfer speed
+    #                        between the data server and the GPU host should be high enough, else the procedure becomes
+    #                        very slow.
     def launch(self, path_data, path_target, objlist_train, objlist_valid):
-        # Build network:
-        self.net.compile(optimizer=self.optimizer, loss=losses.tversky_loss, metrics=['accuracy'])
-
         # Load whole dataset:
         if self.flag_direct_read == False:
             self.display('Loading dataset ...')
@@ -219,14 +240,25 @@ class Train(DeepFinder):
         self.display("Model took %0.2f seconds to train" % np.sum(process_time))
         self.net.save(self.path_out+'net_weights_FINAL.h5')
 
-    # This function generates training batches:
-    #   - Data and target patches are sampled, in order to avoid loading whole tomograms.
+    # Generates batches for training and validation. In this version, the dataset is not loaded into memory. Only the
+    # current batch content is loaded into memory, which is useful when memory is limited.
+    # Is called when self.flag_direct_read=True
+    # !! For now only works for h5 files !!
+    # Batches are generated as follows:
     #   - The positions at which patches are sampled are determined by the coordinates contained in the object list.
     #   - Two data augmentation techniques are applied:
     #       .. To gain invariance to translations, small random shifts are added to the positions.
     #       .. 180 degree rotation around tilt axis (this way missing wedge orientation remains the same).
     #   - Also, bootstrap (i.e. re-sampling) can be applied so that we have an equal amount of each macromolecule in each batch.
-    #     This is usefull when a class is under-represented.
+    #     This is usefull when a class is under-represented. Is applied when self.flag_batch_bootstrap=True
+    # INPUTS:
+    #   path_data: list of strings '/path/to/tomo.h5'
+    #   path_target: list of strings '/path/to/target.h5'
+    #   batch_size: int
+    #   objlist: list of dictionnaries
+    # OUTPUT:
+    #   batch_data: numpy array [batch_idx, z, y, x, channel] in our case only 1 channel
+    #   batch_target: numpy array [batch_idx, z, y, x, class_idx] is one-hot encoded
     def generate_batch_direct_read(self, path_data, path_target, batch_size, objlist=None):
         p_in = np.int(np.floor(self.dim_in / 2))
 
@@ -253,11 +285,11 @@ class Train(DeepFinder):
 
             # Load data and target patches:
             h5file = h5py.File(path_data[tomoID], 'r')
-            patch_data = h5file['dataset'][x - p_in:x + p_in, y - p_in:y + p_in, z - p_in:z + p_in]
+            patch_data = h5file['dataset'][z - p_in:z + p_in, y - p_in:y + p_in, x - p_in:x + p_in]
             h5file.close()
 
             h5file = h5py.File(path_target[tomoID], 'r')
-            patch_target = h5file['dataset'][x - p_in:x + p_in, y - p_in:y + p_in, z - p_in:z + p_in]
+            patch_target = h5file['dataset'][z - p_in:z + p_in, y - p_in:y + p_in, x - p_in:x + p_in]
             h5file.close()
 
             # Process the patches in order to be used by network:
@@ -275,6 +307,17 @@ class Train(DeepFinder):
 
         return batch_data, batch_target
 
+    # Generates batches for training and validation. In this version, the whole dataset has already been loaded into
+    # memory, and batch is sampled from there. Apart from that does the same as above.
+    # Is called when self.flag_direct_read=False
+    # INPUTS:
+    #   data: list of numpy arrays
+    #   target: list of numpy arrays
+    #   batch_size: int
+    #   objlist: list of dictionnaries
+    # OUTPUT:
+    #   batch_data: numpy array [batch_idx, z, y, x, channel] in our case only 1 channel
+    #   batch_target: numpy array [batch_idx, z, y, x, class_idx] is one-hot encoded
     def generate_batch_from_array(self, data, target, batch_size, objlist=None):
         p_in = np.int(np.floor(self.dim_in / 2))
 
@@ -304,8 +347,8 @@ class Train(DeepFinder):
             x, y, z = core_utils.get_patch_position(tomodim, p_in, objlist[index], self.Lrnd)
 
             # Get patch:
-            patch_data = sample_data[x - p_in:x + p_in, y - p_in:y + p_in, z - p_in:z + p_in]
-            patch_batch = sample_target[x - p_in:x + p_in, y - p_in:y + p_in, z - p_in:z + p_in]
+            patch_data  = sample_data[  z - p_in:z + p_in, y - p_in:y + p_in, x - p_in:x + p_in]
+            patch_batch = sample_target[z - p_in:z + p_in, y - p_in:y + p_in, x - p_in:x + p_in]
 
             # Process the patches in order to be used by network:
             patch_data = (patch_data - np.mean(patch_data)) / np.std(patch_data)  # normalize
@@ -344,6 +387,8 @@ class Segment(DeepFinder):
     #   weights_path: path to the .h5 file containing the network weights obtained by the training procedure (string)
     # OUTPUT:
     #   predArray: a numpy array containing the predicted score maps.
+    # Note: in this function x,y,z is actually z,y,x. Has no incidence when used. This note is for someone who wants
+    #       to understand this code.
     def launch(self, dataArray):
         dataArray = (dataArray[:] - np.mean(dataArray[:])) / np.std(dataArray[:])  # normalize
         dataArray = np.pad(dataArray, self.pcrop, mode='constant', constant_values=0)  # zeropad

examples/training/step2_launch_training.py

@@ -27,6 +27,7 @@ Nclass = 13
 # Initialize training task:
 trainer = df.Train(Ncl=Nclass)
 trainer.path_out         = 'out/' # output path
+trainer.h5_dset_name     = 'dataset' # if training data is stored as h5, you can specify the h5 dataset
 trainer.dim_in           = 56 # patch size
 trainer.batch_size       = 25
 trainer.epochs           = 2

pyqt/generate_target/gui_target.py

@@ -95,7 +95,7 @@ class TargetGenerationWindow(QtWidgets.QMainWindow, Ui_MainWindow):
             tbuild.display('Loading initial volume ...')
             vol_initial = utils.read_array(path_initial_vol)
         else:
-            vol_initial = np.zeros((dim_x, dim_y, dim_z))
+            vol_initial = np.zeros((dim_z, dim_y, dim_x))
 
         if strategy == 'Shapes':
             target = tbuild.generate_with_shapes(objl, vol_initial, param_list)

pyqt/segment/gui_segment.ui

@@ -45,7 +45,7 @@
        <item row="1" column="1">
         <widget class="QLineEdit" name="le_path_weights">
          <property name="text">
-          <string>/Users/emoebel/serpico-fs2/github/deep-finder/examples/training/params_model_FINAL.h5</string>
+          <string>/Users/emoebel/serpico-fs2/github/deep-finder/examples/analyze/in/net_weights_FINAL.h5</string>
          </property>
         </widget>
        </item>

utils.py

@@ -136,7 +136,7 @@ def bin_array(array):
 #   orient: list of Euler angles (phi,psi,the) as defined in PyTOM
 # OUTPUT:
 #   arrayR: rotated 3D numpy array
-def rotate_array(array, orient):
+def rotate_array(array, orient): # TODO move to core_utils?
     phi = orient[0]
     psi = orient[1]
     the = orient[2]
@@ -192,7 +192,7 @@ def rotate_array(array, orient):
 #   R  : radius of the sphere (in voxels)
 # OUTPUT:
 #   sphere: 3D numpy array where '1' is 'sphere' and '0' is 'no sphere'
-def create_sphere(dim, R):
+def create_sphere(dim, R): # TODO move to core_utils?
     C = np.floor((dim[0]/2, dim[1]/2, dim[2]/2))
     x,y,z = np.meshgrid(range(dim[0]),range(dim[1]),range(dim[2]))