text: #!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  """
  Created on June 20th 17:28:17 2021
  
  @author: yaleman
  """
  """
  % Syntax :
  % 
  % Description
  %
  %
  % See also: 
  %__________________________________________________
  % Authors: Yasser Aleman Gomez
  % Radiology Department
  % CHUV, Lausanne
  % Created on %(date)s
  % Version $1.0
  
  %% ======================= Importing Libraries ============================== %
  """
  import numpy as np
  import nibabel as nib
  import h5py
  import time
  import os
  from pathlib import Path
  import argparse
  import sys
  import csv
  
  class SmartFormatter(argparse.HelpFormatter):
  
      def _split_lines(self, text, width):
          if text.startswith('R|'):
              return text[2:].splitlines()
          # this is the RawTextHelpFormatter._split_lines
          return argparse.HelpFormatter._split_lines(self, text, width)
  
  
  def _build_args_parser():
      formatter = lambda prog: argparse.HelpFormatter(prog, max_help_position=52)
  
      from argparse import ArgumentParser
  
      p = argparse.ArgumentParser(formatter_class=SmartFormatter, description='\n Help \n')
  
      requiredNamed = p.add_argument_group('Required arguments')
      requiredNamed.add_argument('--method', '-m', action='store', choices=['connectivity', 'consfile'], required=True,
                  metavar='METHOD', type=str, nargs=1, default='connectivity',
                  help="R|Method selection (connectivity or consfile). Another file should be supplied if connectivity is \n"
                       "                                             selected as method. This file should be supplied \n"
                       "                                             through the flags --map -p. \n"
                       "\n"
                       "   connectivity: Compute the mean value of a supplied scalar map along all the bundles included in the \n"
                       "                 atlas. It outputs a connectivity matrix where the connection strength is the mean \n"
                       "                 value of the scalar map along each bundle. \n"
                       " \n"
                       "    consfile   : Extract and save some especific bundles. The selected bundles should be supplied \n"
                       "                 through a nifti-1 file or a csv text file using the flags --bfile -bf. If a scalar \n"
                       "                 map is supplied, a table with the mean values along these bundles is stored as well.\n"
                       "                 - If the file is a Nifti image the bundles intercepting its non-zero values will \n"
                       "                 be extracted and saved. \n"
                       "                 - If the file is a readable text file, each row of the txt file should contain \n"
                       "                 the source and target regions separated by a comma (ie. 10,57). \n"
                       "\n")
      requiredNamed.add_argument('--out', '-o', action='store', required=True, metavar='OUT', type=str, nargs=1,
                  help= "R| Output basename. \n"
                        "                    If connectivity was selected as method then the output filename will be: \n"
                        "                    <outbasename>-connmat-scaleId' \n"
                        "\n"
                        "                    Otherwise, it will be: \n"
                        "                    <outbasename>-<sourceId_targetId>-scaleId')\n")
      requiredNamed.add_argument('--scale', '-s', action='store', required=False, choices=['scale1', 'scale2', 'scale3', 'scale4'],
                 metavar='SCALE', type=str, nargs=1,
                 help= "R| Scale Id (scale1, scale2, scale3 or scale4).\n", default=['scale1'])
      requiredNamed.add_argument('--bfile', '-bf', action='store', required=False, metavar='FILE', type=str, nargs=1,
                  help="R| It can be either a binary mask or a Comma Separated text file.\n"
                       "     - If this file is a Nifti image the bundles intercepting its non-zero values will be \n"
                       "       extracted and saved \n"
                       "     - If the flie is a text file, each row of the txt file should contain the source and target \n"
                       "       regions separated by a comma (ie. 10,57). The first row of this file will be ignored. \n"
                       "\n", default='None')
      requiredNamed.add_argument('--map', '-p', action='store', required=False,
                  metavar='MAP', type=str, nargs=1,
                  help="R| Scalar map to compute its mean value along the selected bundles.\n"
                      "\n" , default='None')
      requiredNamed.add_argument('--subth', '-st', action='store', required=False,
                                 metavar='SUBTHRESH', type=float, nargs=1,
                                 help="R| Subject-level consistency threshold in percentage (0-100). Bundles appearing in \n"
                                      " a percentage of subjects lower than the selected threshold will not be taken into \n"
                                      " account (default = 0). \n"
                                      "\n", default=[0])
      requiredNamed.add_argument('--voxth', '-vt', action='store', required=False,
                                 metavar='VOXTHRESH', type=float, nargs=1,
                                 help="R| Voxel-level consistency threshold. Voxels with probability (0-1) values lower \n"
                                      " than the supplied threshold will be set to 0 in the selected bundles (default = 0).\n"
                                      "\n", default=[0])
      requiredNamed.add_argument('--extract', '-e', action='store_true', required=False,
                                 help="R| Save individual bundles in Nifti-1 format. \n"
                                      "\n")
      requiredNamed.add_argument('--all', '-a', action='store_true', required=False,
                                 help="R| Save individual bundles intercepting any of the regions of interest inside the mask\n"
                                      "\n")
      requiredNamed.add_argument('--collapse', '-c', action='store_true', required=False,
                                 help="R| Collapse the bundles into a 4D color-coded Nifti-1 file. \n"
                                      "\n")
      requiredNamed.add_argument('--force', '-f', action='store_true', required=False,
                                 help="R| Overwrite the results. \n"
                                      "\n")
      p.add_argument('--verbose', '-v', action='store', required=False,
                  type=int, nargs=1,
                  help='verbosity level: 1=low; 2=debug')
  
      args   = p.parse_args()
      method = args.method[0]
      map    = args.map[0]
      file   = args.bfile[0]
  
      # boolmap = bool(map) + bool(event_id) + bool(broker_id) + bool(product_id) + bool(event_keyword) + bool(metadata)
      if method == 'connectivity' and not os.path.isfile(map):
          print("\n")
          print("Please, supply a valid scalar map using the flags --map or -p.")
          p.print_help()
          sys.exit()
      elif method != 'connectivity' and not os.path.isfile(file):
          print(" ")
          print(" ")
          print("Please, supply a valid txt or binary mask to select the bundles. Use the flags --bfile or -f.")
          print(" ")
          print(" ")
  
          p.print_help()
          sys.exit()
  
      return p
  
  def save_matrix(connMat, matFilename, lnames, regCodes, regCoords, regColors):
      """
      Save the connectivity matrix to a HDF5 file.
      @params:
          connMat        - Required  : Numpy array with the connectivity matrix (Float)
          matFilename    - Required  : Output filename for the connectivity matrix (Str)
          lnames         - Required  : Region names (Str)
          regCodes       - Required  : Region codes (position in the matrix) (Int)
          regCoords      - Required  : Coordinates of the regions center-of mass in mm (MNI space)(Float)
          regColors      - Required  : Region RGB (red, green and blue) colors triplets. It can be used for networks visualization.(Float)
      """
  
      # import networkx as nx
      import h5py
      # G = nx.from_numpy_matrix(connMat)
      # nx.write_gpickle(G, matFilename)
  
      indnan = np.isnan(connMat)
      connMat[indnan] = 0
  
      hm = h5py.File(matFilename, 'w')
      # Saving the Connectivity Matrix
      g1 = hm.create_group('connmat')
      g1.create_dataset('matrix', data=connMat)
      g1.create_dataset('gmregions', data=np.array(lnames, dtype='S'))
      g1.create_dataset('gmcodes', data=regCodes)
      g1.create_dataset('gmcoords', data=regCoords)
      g1.create_dataset('gmcolors', data=regColors)
      hm.close()
  
      return matFilename
  
  # Print iterations progress
  def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1, length = 100, fill = '█', printEnd = "\r"):
      """
      Call in a loop to create terminal progress bar
      @params:
          iteration   - Required  : current iteration (Int)
          total       - Required  : total iterations (Int)
          prefix      - Optional  : prefix string (Str)
          suffix      - Optional  : suffix string (Str)
          decimals    - Optional  : positive number of decimals in percent complete (Int)
          length      - Optional  : character length of bar (Int)
          fill        - Optional  : bar fill character (Str)
          printEnd    - Optional  : end character (e.g. "\r", "\r\n") (Str)
      """
      percent = ("{0:." + str(decimals) + "f}").format(100 * (iteration / float(total)))
      filledLength = int(length * iteration // total)
      bar = fill * filledLength + '-' * (length - filledLength)
      print(f'\r{prefix} |{bar}| {percent}% {suffix}', end = printEnd)
      # Print New Line on Complete
      if iteration == total:
          print()
  
  def bundles_from_mask(hfName, maskFilename, outBasename, scaleId, mapFilename, voxth=0, subth=0, boolall=False, boolbund=False, boolcoll=False):
      """
      Compute mean, median and a standard deviation of a certain scalar map along the bundles reaching specified ROIs
      @params:
          hfName            - Required  : h5 filename (Str)
          maskFilename      - Required  : Nifti-1 filename containing the ROIs. Each ROI should have a different value (Str)
          outBasename       - Required  : Output basename (Str)
          scaleId           - Required  : ScaleId: scale1, scale2, scale3 or scale4  (Str)
          mapFilename       - Required  : Scalar map filename (Int)
          voxth             - Optional  : Voxel-wise consistency threshold (Float)
          subth             - Optional  : Subject-wise consistency threshold (Float)
          boolall           - Optional  : Boolean variable to save all the bundles reaching any of the selected ROIs(Bool)
          boolbund          - Optional  : Boolean variable to save the selected bundles (Bool)
          boolcoll          - Optional  : Boolean variable to collapse the selected bundles (Bool)
      """
  
      outDir = os.path.dirname(outBasename)
      fname  = os.path.basename(outBasename)
  
      # Creating the folder if it does not exist
      if not os.path.isdir(outDir):
          path = Path(outDir)
          try:
              path.mkdir()
          except OSError:
              print("Failed to make nested directory")
          else:
              print("Nested directory made")
  
      # Reading the mask image
      maskI = nib.load(maskFilename)
      dataP = maskI.get_fdata()
  
      if boolall:
          dataP[dataP != 0] = 1
  
      nrois_mask = np.unique(dataP)
      nrois_mask = nrois_mask[nrois_mask != 0]
  
      # Reading the HDF5 file
      hf        = h5py.File(hfName, 'r')
      tempVar   = hf.get('header/affine')
      affine    = np.array(tempVar)
      tempVar   = hf.get('header/dim')
      dim       = np.array(tempVar)
      tempVar   = hf.get('header/gmregions')
      stnames   = np.array(tempVar)
      tempVar   = hf.get('header/nsubjects')
      nsubj     = np.array(tempVar)
      tempVar   = hf.get('header/gmcolors')
      gmcolors  = np.array(tempVar)
  
      # Reading the consistency matrix. It can be used for thresholding
      tempVar   = hf.get('matrices/consistency')
      consistM  = np.array(tempVar)/nsubj
  
      # tempVar = hf.get('header/scale')
      # scaleId = np.array(tempVar)
      if boolcoll:
          collMask  = np.zeros((dim[0], dim[1], dim[2]), dtype='float')  # Collapsing all the bundles
          collMaskc = np.zeros((dim[0], dim[1], dim[2], 3), dtype='float')
  
  
  
      # Bundles to be used
      bundloc   = np.argwhere(np.triu(consistM, 1) >= subth/100)
      X         = bundloc[:, 0] + 1
      Y         = bundloc[:, 1] + 1
  
  
      # Reading the scalar map
      booltab     = False
      if os.path.isfile(mapFilename):
          mapI    = nib.load(mapFilename)
          dataM   = mapI.get_fdata()
          booltab = True
  
      tabList   = [['BundleId', 'SourceROI', 'TargetROI', 'Consistency' , 'MeanValue', 'StdValue', 'MedianValue']]
      tic       = time.perf_counter()
  
      # Reading the bundles keys
      # tempVar       = hf.get('atlas')
      # bundIds       = list(tempVar.keys())
      # tbund         = len(bundIds)
      tbund     = len(X)
      # Loop along the bundles
      for i in np.arange(0, tbund - 1):
          # tempVar = hf.get("atlas/" + bundIds[i])
          tempS      = str(X[i]) + "_" + str(Y[i])
          tempVar    = hf.get("atlas/" + tempS)
          voxIds     = np.array(tempVar)
  
          bdcolor    = np.mean(gmcolors[bundloc[i, :],:],axis=0)/255
  
          # Creating the empty volume
          array_data = np.zeros((dim[0], dim[1], dim[2]), dtype=float)
  
          # Assigning the probability
          array_data[voxIds[:, 0], voxIds[:, 1], voxIds[:, 2]] = voxIds[:, 3] / nsubj
  
          # Thresholding using the concurrency value
          if voxth > 0:
              array_data[array_data < voxth] = 0
  
          # Interception between the mask and the probabilistic atlas
          tempM      = np.multiply(array_data != 0, dataP)
  
          nrois_bund = np.unique(tempM)
          nrois_bund = nrois_bund[nrois_bund != 0]
  
          savebund = False
          if len(nrois_mask) == 1 and len(nrois_bund) == 1:
              savebund = True
          elif len(nrois_mask) != 1 and len(nrois_bund) > 1:
              savebund = True
  
          if savebund:
              if boolcoll:
                  collMask            = collMask  + np.asarray(array_data >0)
                  collMaskc[:,:,:,0]  = collMaskc[:,:,:,0] + array_data*bdcolor[0]
                  collMaskc[:,:,:,1]  = collMaskc[:,:,:,1] + array_data*bdcolor[1]
                  collMaskc[:,:,:,2]  = collMaskc[:,:,:,2] + array_data*bdcolor[2]
  
              # Saving the bundle image
              if boolbund:
                  bundFilename   = os.path.join(outDir, fname + "-fromMask-" + tempS + "-" + scaleId + ".nii.gz")
                  array_img      = nib.Nifti1Image(array_data, affine)
                  nib.save(array_img, bundFilename)
  
              # Creating the table
              if booltab:
                  # tempS = bundIds[i]
                  # roiIds = tempS.split('_')
  
                  mVal    = np.mean(dataM[array_data >= voxth])
                  sVal    = np.std(dataM[array_data >= voxth])
                  dVal    = np.median(dataM[array_data >= voxth])
                  sourROI = stnames[int(X[i])-1].decode('UTF-8')
                  targROI = stnames[int(Y[i])-1].decode('UTF-8')
                  cVal    = consistM[X[i]-1, Y[i]-1]*100
                  tabList.append([tempS, sourROI, targROI, cVal, mVal, sVal, dVal])
  
          printProgressBar(i + 1, tbund , 'Processing bundle: ' + '(' + str(i+1) + '/' + str(tbund) + ')')
  
      # Saving the table
      if booltab:
          tabFilename     = outBasename + "-table-fromMask-" + scaleId + ".csv"
          with open(tabFilename, mode='w') as bundTab:
              bund = csv.writer(bundTab, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
              for i in np.arange(len(tabList)):
                  bund.writerow(tabList[i])
      # Saving the collapsed image
      if boolcoll:
          collFilename           = os.path.join(outDir, fname + "-fromMask-bundcollap-" + scaleId + ".nii.gz")
          # collMask        = collMask != 0
          # array_img              = nib.Nifti1Image(collMask.astype(int), affine)
          # nib.save(array_img, collFilename)
          collMaskc[:, :, :, 0] = np.divide(collMaskc[:, :, :, 0], collMask, where=True)
          collMaskc[:, :, :, 1] = np.divide(collMaskc[:, :, :, 1], collMask, where=True)
          collMaskc[:, :, :, 2] = np.divide(collMaskc[:, :, :, 2], collMask, where=True)
          collMaskc[np.isnan(collMaskc)] = 0
  
          imgcoll               = nib.Nifti1Image(collMaskc, affine)
          nib.save(imgcoll, collFilename)
  
      # Timer
      toc                 = time.perf_counter()
      print(f"Processed in {toc - tic:0.4f} seconds")
  
  def bundles_from_textfile(hfName, textFilename, outBasename, scaleId, mapFilename, voxth=0, subth=0, boolbund=False, boolcoll=False):
      """
      Compute mean, median and a standard deviation of a certain scalar map along the bundles specified in a csv
      file.
      @params:
          hfName            - Required  : h5 filename (Str)
          textFilename      - Required  : Comma Separated Value file with the bundles of interest (i.e. 10_57, bundle between regions 10 and 57) (Str)
          outBasename       - Required  : Output basename (Str)
          scaleId           - Required  : ScaleId: scale1, scale2, scale3 or scale4  (Str)
          mapFilename       - Required  : Scalar map filename (Int)
          voxth             - Optional  : Voxel-wise consistency threshold (range: 0-1)(Float)
          subth             - Optional  : Subject-wise consistency threshold in precentage (range: 0-100)(Float)
          boolbund          - Optional  : Boolean variable to save the selected bundles (Bool)
          boolcoll          - Optional  : Boolean variable to collapse the selected bundles (Bool)
      """
  
      outDir = os.path.dirname(outBasename)
      fname = os.path.basename(outBasename)
  
      # Creating the folder if it does not exist
      if not os.path.isdir(outDir):
          path = Path(outDir)
          try:
              path.mkdir()
          except OSError:
              print("Failed to make nested directory")
          else:
              print("Nested directory made")
  
      # Reading the txt file
      bundIds = np.loadtxt(textFilename, skiprows=1, delimiter=',', dtype='u4', ndmin=2)
      bundIds = np.sort(bundIds, axis=1)
  
      # Reading the HDF5 file
      hf = h5py.File(hfName, 'r')
      tempVar = hf.get('header/affine')
      affine = np.array(tempVar)
      tempVar = hf.get('header/dim')
      dim = np.array(tempVar)
      tempVar = hf.get('header/gmregions')
      stnames = np.array(tempVar)
      tempVar = hf.get('header/nsubjects')
      nsubj = np.array(tempVar)
      tempVar   = hf.get('header/gmcolors')
      gmcolors  = np.array(tempVar)
  
      # Reading the consistency matrix. It can be used for thresholding
      tempVar  = hf.get('matrices/consistency')
      consistM = np.array(tempVar) / nsubj
  
      if boolcoll:
          collMask  = np.zeros((dim[0], dim[1], dim[2]), dtype='float')  # Collapsing all the bundles
          collMaskc = np.zeros((dim[0], dim[1], dim[2], 3), dtype='float')
  
  
      # If the csv file contains only one column then all the bundles involving the selected regions are extracted
      if np.shape(bundIds)[1] == 1:
          nstr     = len(consistM)
          sbundIds = np.zeros([1, 2], dtype='u4')
          for i in np.arange(0, len(bundIds)):
              a1 = np.ones([nstr, 1], dtype='u4') * bundIds[i]
              a2 = np.linspace(1, nstr, nstr, dtype='u4')
              a2 = a2[:, None]
              out_mat = np.concatenate((a1, a2), axis=1)
              sbundIds = np.concatenate((sbundIds, out_mat), axis=0)
          sbundIds = np.sort(sbundIds, axis=1)
          t = sbundIds[:, 0] - sbundIds[:, 1] != 0
          bundIds = sbundIds[t, :]
  
      # Bundles to be used
      bundloc  = np.argwhere(np.triu(consistM, 1) >= subth / 100)
      X        = bundloc[:, 0] + 1
      Y        = bundloc[:, 1] + 1
  
  
      # Reading the scalar map
      booltab = False
      if os.path.isfile(mapFilename):
          mapI = nib.load(mapFilename)
          dataM = mapI.get_fdata()
          booltab = True
  
      # Loop along the bundles
      # p = ProgressBar(len(bundIds))
      tabList   = [['BundleId', 'SourceROI', 'TargetROI', 'Consistency' ,'MeanValue', 'StdValue', 'MedianValue']]
      tic = time.perf_counter()
      tbund = len(bundIds)
  
      for i in np.arange(tbund):
          t = gmcolors[bundIds[i, :]-1, :]
          bdcolor = np.mean(t, axis=0) / 255
  
          # Reading the selected bundles
          tempVar = hf.get("atlas/" + str(bundIds[i, 0]) + "_" + str(bundIds[i, 1]))
  
          if tempVar != None:
              voxIds = np.array(tempVar)
  
              # Creating the empty volume
              array_data = np.zeros((dim[0], dim[1], dim[2]), dtype=float)
  
              # Assigning the probability
              array_data[voxIds[:, 0], voxIds[:, 1], voxIds[:, 2]] = voxIds[:, 3] / nsubj
  
              # Thresholding using the currency value
              if voxth > 0:
                  array_data[array_data < voxth] = 0
  
              # Writing the nifti file
              if np.any(array_data):
                  collMask            = collMask  + np.asarray(array_data >0)
                  collMaskc[:,:,:,0]  = collMaskc[:,:,:,0] + array_data*bdcolor[0]
                  collMaskc[:,:,:,1]  = collMaskc[:,:,:,1] + array_data*bdcolor[1]
                  collMaskc[:,:,:,2]  = collMaskc[:,:,:,2] + array_data*bdcolor[2]
                  if boolbund:
                      # Saving only the bundles intercepting the mask
                      bundFilename = os.path.join(outDir, fname + "-fromCSV-" + str(bundIds[i, 0]) + "_" + str(bundIds[i, 1]) + "-" + scaleId + ".nii.gz")
                      array_img = nib.Nifti1Image(array_data, affine)
                      nib.save(array_img, bundFilename)
  
                  # Creating the table
                  if booltab:
                      mVal      = np.mean(dataM[array_data >= voxth])
                      sVal      = np.std(dataM[array_data >= voxth])
                      dVal      = np.median(dataM[array_data >= voxth])
                      sourROI   = stnames[bundIds[i, 0]-1].decode('UTF-8')
                      targROI   = stnames[bundIds[i, 1]-1].decode('UTF-8')
                      cVal = consistM[bundIds[i, 0] - 1, bundIds[i, 1]  - 1] * 100
                      tabList.append([str(bundIds[i, 0]) + "_" + str(bundIds[i, 1]), sourROI, targROI, cVal, mVal, sVal, dVal])
  
              printProgressBar(i+1, tbund, 'Processing bundle: ' + '(' + str(i + 1) + '/' + str(tbund) + ')')
  
      # Saving the table
      if booltab:
          tabFilename = outBasename + "-table-fromCSV-" + scaleId + ".csv"
          with open(tabFilename, mode='w') as bundTab:
              bund = csv.writer(bundTab, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
              for i in np.arange(len(tabList)):
                  bund.writerow(tabList[i])
  
      # Saving the collapsed image
      if boolcoll:
          collFilename = os.path.join(outDir, fname + "-fromCSV-bundcollap-" + scaleId + ".nii.gz")
          collMaskc[:, :, :, 0] = np.divide(collMaskc[:, :, :, 0], collMask, where=True)
          collMaskc[:, :, :, 1] = np.divide(collMaskc[:, :, :, 1], collMask, where=True)
          collMaskc[:, :, :, 2] = np.divide(collMaskc[:, :, :, 2], collMask, where=True)
          collMaskc[np.isnan(collMaskc)] = 0
  
          imgcoll               = nib.Nifti1Image(collMaskc, affine)
          nib.save(imgcoll, collFilename)
  
      # Timer
      toc = time.perf_counter()
      print(f"Processed in {toc - tic:0.4f} seconds")
  
  def bundles_to_connectivity(hfName, outBasename, scaleId, mapFilename, voxth=0, subth=0, boolbund=False, boolforce=False):
      """
       Compute connectivity matrix for a specified scale where the streng of connection is mean, median and standard deviation of a certain scalar map .
       @params:
           hfName            - Required  : h5 filename (Str)
           textFilename      - Required  : Comma Separated Value file with the bundles of interest (i.e. 10_57, bundle between regions 10 and 57) (Str)
           outBasename       - Required  : Output basename (Str)
           scaleId           - Required  : ScaleId: scale1, scale2, scale3 or scale4  (Str)
           mapFilename       - Required  : Scalar map filename (Int)
           voxth             - Optional  : Voxel-wise consistency threshold (range: 0-1)(Float)
           subth             - Optional  : Subject-wise consistency threshold in precentage (range: 0-100)(Float)
           boolbund          - Optional  : Boolean variable to save the selected bundles (Bool)
           boolforce         - Optional  : Rewrite the connectivity matrices (Bool)
       """
      outDir = os.path.dirname(outBasename)
      fname = os.path.basename(outBasename)
  
      # Creating the folder if it does not exist
      if not os.path.isdir(outDir):
          path = Path(outDir)
          try:
              path.mkdir()
          except OSError:
              print("Failed to make nested directory")
          else:
              print("Nested directory made")
  
      outBasename + "-connmat-mean-" + scaleId + ".h5"
  
  
      if os.path.isfile(outBasename + "-connmat-mean-" + scaleId + ".h5") and boolforce == True:
          print(" The file " + outBasename + "-connmat-mean-" + scaleId + ".h5" + " will be overwritten.")
          os.remove(outBasename + "-connmat-mean-" + scaleId + ".h5")
          print("  ")
      if os.path.isfile(outBasename + "-connmat-std-" + scaleId + ".h5") and boolforce == True:
          print(" The file " + outBasename + "-connmat-std-" + scaleId + ".h5" + " will be overwritten.")
          os.remove(outBasename + "-connmat-std-" + scaleId + ".h5")
          print("  ")
      if os.path.isfile(outBasename + "-connmat-median-" + scaleId + ".h5") and boolforce == True:
          print(" The file " + outBasename + "-connmat-median-" + scaleId + ".h5" + " will be overwritten.")
          os.remove(outBasename + "-connmat-median-" + scaleId + ".h5")
          print("  ")
  
      # Reading the HDF5 file
      hf        = h5py.File(hfName, 'r')
      tempVar   = hf.get('header/affine')
      affine    = np.array(tempVar)
      tempVar   = hf.get('header/dim')
      dim       = np.array(tempVar)
  
      # Reading the fields to plot the network
      tempVar   = hf.get('header/gmregions')
      stnames   = np.array(tempVar)
      tempVar = hf.get('header/gmcodes')
      regCodes = np.array(tempVar)
      tempVar = hf.get('header/gmcoords')
      regCoords = np.array(tempVar)
      tempVar = hf.get('header/gmcolors')
      regColors = np.array(tempVar)
  
      # Reading the number of subjects used to create the atlas
      tempVar   = hf.get('header/nsubjects')
      nsubj     = np.array(tempVar)
  
      # Reading the consistency matrix. It can be used for thresholding
      tempVar   = hf.get('matrices/consistency')
      consistM  = np.array(tempVar)/nsubj
  
      # Bundles to be used
      bundloc = np.argwhere(np.triu(consistM, 1) > subth/100)
      X = bundloc[:, 0] + 1
      Y = bundloc[:, 1] + 1
  
      # tempVar = hf.get('header/scale')
      # scaleId = np.array(tempVar)
  
      # Reading the scalar map
      booltab = False
      if os.path.isfile(mapFilename):
          mapI = nib.load(mapFilename)
          dataM = mapI.get_fdata()
          booltab = True
  
      # Creating empty connectivity matrices
      nstr        = len(stnames)
      meanMat     = np.zeros([nstr, nstr])  # Number of fibers
      stdMat      = np.zeros([nstr, nstr])  # Fiber Length
      medMat      = np.zeros([nstr, nstr])  # Number of Points. Important to compute the mean along fibers
  
  
      tic     = time.perf_counter()
  
      # Reading the bundles keys
      # tempVar       = hf.get('atlas')
      # bundIds       = list(tempVar.keys())
      # tbund         = len(bundIds)
  
      tbund = len(X)
      # Loop along the bundles
      for i in np.arange(0, tbund - 1):
          # tempS   = bundIds[i]
          tempS   = str(X[i]) + "_" + str(Y[i])
          tempVar = hf.get("atlas/" + tempS)
          voxIds  = np.array(tempVar)
          # roiIds  = tempS.split('_')
  
          # Creating the empty volume
          array_data = np.zeros((dim[0], dim[1], dim[2]), dtype=float)
  
          # Assigning the probability
          array_data[voxIds[:, 0], voxIds[:, 1], voxIds[:, 2]] = voxIds[:, 3] / nsubj
  
          # # Thresholding using the currency value
          # if voxth > 0:
          #     array_data[array_data < voxth] = 0
  
          # Writing the nifti file
          if np.any(array_data):
  
              if boolbund:
                  # Saving only the bundles intercepting the mask
                  # bundFilename = os.path.join(outDir, fname + "-atlas-" + str(roiIds[0]) + "_" + str(roiIds[1]) + "-" + scaleId + ".nii.gz")
                  bundFilename = os.path.join(outDir, fname + "-atlas-" + str(X[i]) + "_" + str(Y[i]) + "-" + scaleId + ".nii.gz")
                  array_img = nib.Nifti1Image(array_data, affine)
                  nib.save(array_img, bundFilename)
  
              # Creating the connectivity matrices
              # meanMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1] = np.mean(dataM[array_data >= voxth])
              # meanMat[int(roiIds[1]) - 1, int(roiIds[0]) - 1] = meanMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1]
              meanMat[int(X[i]) - 1, int(Y[i]) - 1] = np.mean(dataM[array_data >= voxth])
              meanMat[int(Y[i]) - 1, int(X[i]) - 1] = meanMat[int(X[i]) - 1, int(Y[i]) - 1]
  
              # stdMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1] = np.std(dataM[array_data >= voxth])
              # stdMat[int(roiIds[1]) - 1, int(roiIds[0]) - 1] = stdMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1]
              stdMat[int(X[i]) - 1, int(Y[i]) - 1] = np.mean(dataM[array_data >= voxth])
              stdMat[int(Y[i]) - 1, int(X[i]) - 1] = stdMat[int(X[i]) - 1, int(Y[i]) - 1]
  
              # medMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1] = np.median(dataM[array_data >= voxth])
              # medMat[int(roiIds[1]) - 1, int(roiIds[0]) - 1] = medMat[int(roiIds[0]) - 1, int(roiIds[1]) - 1]
              medMat[int(X[i]) - 1, int(Y[i]) - 1] = np.mean(dataM[array_data >= voxth])
              medMat[int(Y[i]) - 1, int(X[i]) - 1] = medMat[int(X[i]) - 1, int(Y[i]) - 1]
  
          printProgressBar(i+1, tbund, 'Processing bundle: ' + '(' + str(i+1) + '/' + str(tbund) + ')')
  
      # Saving the connectivity matrices
  
      # Mean
      meanFilename = outBasename + "-connmat-mean-" + scaleId + ".h5"
      # save_matrix(meanMat, meanFilename, gmcoords, , stnames)
      save_matrix(meanMat, meanFilename, stnames, regCodes, regCoords, regColors)
  
      # Std
      stdFilename = outBasename + "-connmat-std-" + scaleId + ".h5"
      # save_matrix(stdMat, stdFilename)
      save_matrix(stdMat, stdFilename, stnames, regCodes, regCoords, regColors)
  
      # Median
      medFilename = outBasename + "-connmat-median-" + scaleId + ".h5"
      # save_matrix(medMat, medFilename)
      save_matrix(medMat, medFilename, stnames, regCodes, regCoords, regColors)
  
      # Timer
      toc = time.perf_counter()
      print(f"Processed in {toc - tic:0.4f} seconds")
  
  def main():
      import sys
      # 0. Handle inputs
      parser = _build_args_parser()
      args = parser.parse_args()
          
      print(args)
      if args.verbose is not None:
          v = np.int(args.verbose[0])
      else:
          v = 0
          print('- Verbose set to 0\n')
      if v:
          print('\nInputs\n')
  #
  
      # Getting the path of the current running python file
      cwd = os.getcwd()
  
      scaleid   = args.scale[0]
      file      = args.bfile[0]
      fmeth     = args.method[0]
      fout      = args.out[0]
      fmap      = args.map[0]
      voxth     = args.voxth[0]
      subth     = args.subth[0]
      boolall   =  args.all
      boolbund  = args.extract
      boolcoll  = args.collapse
      boolforce = args.force
  
      hfName = os.path.join(cwd,'probconnatlas', "wm.connatlas." + scaleid + ".h5")
  
      if fmeth == 'connectivity':
          bundles_to_connectivity(hfName, fout, scaleid, fmap, voxth, subth, boolbund, boolforce)
  
      elif fmeth == 'consfile':
          # fid = open(file, "r+")
          imask = False
          if file[-7:]== '.nii.gz' or file[-4:] == '.nii':
              imask = True
  
          if imask == True:
              bundles_from_mask(hfName, file, fout, scaleid, fmap, voxth, subth, boolall, boolbund, boolcoll)
          else:
              bundles_from_textfile(hfName, file, fout, scaleid, fmap, voxth, subth, boolbund, boolcoll)
      else:
          print("Please enter a correct method (connectivity or consfile). ")
          sys.exit()
  
  
  if __name__ == "__main__":
      main()
  