diff --git a/environment.yml b/environment.yml index 54d4265..e014b80 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ channels: - conda-forge - bioconda dependencies: - - python=3.7 + - python=3.8 - cython==0.29.24 - pillow==8.3.2 - pypubsub==4.0.3 @@ -17,13 +17,21 @@ dependencies: - scipy==1.7.1 - vtk==9.0.3 - wxpython==4.1.1 + - mido==1.2.10 + - nest-asyncio==1.5.1 - pip - pip: - - python-gdcm==3.0.9.1 + - aioconsole==0.3.2 + - opencv-python==4.5.3.56 - plaidml-keras==0.7.0 - - theano==1.0.5 - - pyacvd==0.2.3 + - python-gdcm==3.0.9.1 + - python-rtmidi==1.4.9 + - python-socketio[client]==5.3.0 + - pyacvd==0.2.7 - pyclaron - polhemusFT - polhemus - pypolaris + - theano==1.0.5 + - uvicorn[standard]==0.15.0 + - https://github.com/baranaydogan/trekker/raw/master/binaries/Trekker-0.9-cp38-cp38-win_amd64.whl diff --git a/invesalius/constants.py b/invesalius/constants.py index 77630cd..3f81dc0 100644 --- a/invesalius/constants.py +++ b/invesalius/constants.py @@ -823,23 +823,35 @@ COIL_ANGLE_ARROW_PROJECTION_THRESHOLD = 5 CAM_MODE = True # Tractography visualization -N_TRACTS = 100 -PEEL_DEPTH = 5 -MAX_PEEL_DEPTH = 30 -SEED_OFFSET = 15 +N_TRACTS = 200 +PEEL_DEPTH = 10 +MAX_PEEL_DEPTH = 40 +SEED_OFFSET = 30 SEED_RADIUS = 1.5 # Increased the default sleep parameter from 0.1 to 0.15 to decrease CPU load during navigation. -SLEEP_NAVIGATION = 0.15 +SLEEP_NAVIGATION = 0.2 SLEEP_COORDINATES = 0.05 SLEEP_ROBOT = 0.01 -BRAIN_OPACITY = 0.5 +BRAIN_OPACITY = 0.6 N_CPU = psutil.cpu_count() -TREKKER_CONFIG = {'seed_max': 1, 'step_size': 0.1, 'min_fod': 0.1, 'probe_quality': 3, - 'max_interval': 1, 'min_radius_curv': 0.8, 'probe_length': 0.4, - 'write_interval': 50, 'numb_threads': '', 'max_lenth': 200, - 'min_lenth': 20, 'max_sampling_step': 100} +# the max_sampling_step can be set to something different as well. Above 100 is probably not necessary +TREKKER_CONFIG = {'seed_max': 1, + 'step_size': 0.03125, + 'min_fod': 0.05, + 'probe_quality': 3, + 'max_interval': 1, + 'min_radius_curvature': 0.625, + 'probe_length': 0.15625, + 'write_interval': 50, + 'numb_threads': '', + 'max_length': 250, + 'min_length': 10, + 'max_sampling_step': 100, + 'data_support_exponent': 0.5, + 'use_best_init': True, + 'init_max_est_trials': 100} MARKER_FILE_MAGICK_STRING = "##INVESALIUS3_MARKER_FILE_" CURRENT_MARKER_FILE_VERSION = 0 diff --git a/invesalius/control.py b/invesalius/control.py index 5a42e69..eb0d595 100644 --- a/invesalius/control.py +++ b/invesalius/control.py @@ -69,7 +69,7 @@ class Controller(): #DICOM = 1 #TIFF uCT = 2 self.img_type = 0 - self.affine = None + self.affine = np.identity(4) #Init session session = ses.Session() @@ -340,7 +340,7 @@ class Controller(): if proj.affine: self.Slice.affine = np.asarray(proj.affine).reshape(4, 4) else: - self.Slice.affine = None + self.Slice.affine = np.identity(4) Publisher.sendMessage('Update threshold limits list', threshold_range=proj.threshold_range) @@ -623,6 +623,8 @@ class Controller(): Publisher.sendMessage(('Set scroll position', 'SAGITAL'),index=proj.matrix_shape[1]/2) Publisher.sendMessage(('Set scroll position', 'CORONAL'),index=proj.matrix_shape[2]/2) + # TODO: Check that this is needed with the new way of using affine + # now the affine should be at least the identity(4) and never None if self.Slice.affine is not None: Publisher.sendMessage('Enable Go-to-Coord', status=True) else: @@ -731,6 +733,8 @@ class Controller(): proj.level = self.Slice.window_level proj.threshold_range = int(matrix.min()), int(matrix.max()) proj.spacing = self.Slice.spacing + # TODO: Check that this is needed with the new way of using affine + # now the affine should be at least the identity(4) and never None if self.Slice.affine is not None: proj.affine = self.Slice.affine.tolist() diff --git a/invesalius/data/coordinates.py b/invesalius/data/coordinates.py index 7f9aa20..d952929 100644 --- a/invesalius/data/coordinates.py +++ b/invesalius/data/coordinates.py @@ -371,7 +371,7 @@ def DebugCoordRandom(trk_init, trck_id, ref_mode): # # else: - dx = [-70, 70] + dx = [-30, 30] dt = [-180, 180] coord1 = np.array([uniform(*dx), uniform(*dx), uniform(*dx), diff --git a/invesalius/data/coregistration.py b/invesalius/data/coregistration.py index b9948d0..14625c9 100644 --- a/invesalius/data/coregistration.py +++ b/invesalius/data/coregistration.py @@ -427,7 +427,7 @@ class CoordinateCorregistrateNoObject(threading.Thread): # # seed_aux = pos_world.reshape([1, 4])[0, :3] # # seed = seed_aux[np.newaxis, :] # # -# # self.tracts = dtr.compute_tracts(tracker, seed, affine_vtk, True) +# # self.tracts = dtr.compute_and_visualize_tracts(tracker, seed, affine_vtk, True) # # # wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) # wx.CallAfter(Publisher.sendMessage, 'Update cross position', arg=m_img, position=coord) @@ -654,7 +654,7 @@ class CoordinateCorregistrateNoObject(threading.Thread): # # seed_aux = pos_world.reshape([1, 4])[0, :3] # # seed = seed_aux[np.newaxis, :] # -# # self.tracts = dtr.compute_tracts(tracker, seed, affine_vtk, True) +# # self.tracts = dtr.compute_and_visualize_tracts(tracker, seed, affine_vtk, True) # # # wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) # wx.CallAfter(Publisher.sendMessage, 'Update cross position', arg=m_img, position=coord) diff --git a/invesalius/data/styles.py b/invesalius/data/styles.py index 01f7408..89dbab4 100644 --- a/invesalius/data/styles.py +++ b/invesalius/data/styles.py @@ -544,7 +544,7 @@ class TractsInteractorStyle(CrossInteractorStyle): def OnTractsMouseClick(self, obj, evt): # print("Single mouse click") - # self.tracts = dtr.compute_tracts(self.tracker, self.seed, self.left_pressed) + # self.tracts = dtr.compute_and_visualize_tracts(self.tracker, self.seed, self.left_pressed) self.ChangeTracts(True) def OnTractsReleaseLeftButton(self, obj, evt): @@ -554,7 +554,7 @@ class TractsInteractorStyle(CrossInteractorStyle): def ChangeTracts(self, pressed): # print("Trying to compute tracts") - self.tracts = dtr.compute_tracts(self.tracker, self.seed, self.affine_vtk, pressed) + self.tracts = dtr.compute_and_visualize_tracts(self.tracker, self.seed, self.affine_vtk, pressed) # mouse_x, mouse_y = iren.GetEventPosition() # wx, wy, wz = self.viewer.get_coordinate_cursor(mouse_x, mouse_y, self.picker) # px, py = self.viewer.get_slice_pixel_coord_by_world_pos(wx, wy, wz) diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index 955f5b8..a80ea70 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -365,9 +365,16 @@ class SurfaceManager(): if self.convert2inv: # convert between invesalius and world space with shift in the Y coordinate affine = sl.Slice().affine + # TODO: Check that this is needed with the new way of using affine + # now the affine should be at least the identity(4) and never None if affine is not None: - affine[1, -1] -= sl.Slice().spacing[1] * (sl.Slice().matrix.shape[1] - 1) + matrix_shape = sl.Slice().matrix.shape + spacing = sl.Slice().spacing + img_shift = spacing[1] * (matrix_shape[1] - 1) + affine = sl.Slice().affine.copy() + affine[1, -1] -= img_shift affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(affine) + actor.SetUserMatrix(affine_vtk) if overwrite: diff --git a/invesalius/data/tractography.py b/invesalius/data/tractography.py index 950db16..7eba8e1 100644 --- a/invesalius/data/tractography.py +++ b/invesalius/data/tractography.py @@ -29,14 +29,11 @@ import time import numpy as np import queue from invesalius.pubsub import pub as Publisher -from scipy.stats import norm import vtk import invesalius.constants as const import invesalius.data.imagedata_utils as img_utils -import invesalius.project as prj - # Nice print for arrays # np.set_printoptions(precision=2) # np.set_printoptions(suppress=True) @@ -104,75 +101,57 @@ def compute_tubes(trk, direction): return trk_tube -def combine_tracts_root(out_list, root, n_block): +def create_branch(out_list, n_block): """Adds a set of tracts to given position in a given vtkMultiBlockDataSet :param out_list: List of vtkTubeFilters representing the tracts :type out_list: list - :param root: A collection of tracts as a vtkMultiBlockDataSet - :type root: vtkMultiBlockDataSet :param n_block: The location in the given vtkMultiBlockDataSet to insert the new tracts :type n_block: int - :return: The updated collection of tracts as a vtkMultiBlockDataSet - :rtype: vtkMultiBlockDataSet - """ - - # create tracts only when at least one was computed - # print("Len outlist in root: ", len(out_list)) - if not out_list.count(None) == len(out_list): - for n, tube in enumerate(out_list): - root.SetBlock(n_block + n, tube.GetOutput()) - - return root - - -def combine_tracts_branch(out_list): - """Combines a set of tracts in vtkMultiBlockDataSet - - :param out_list: List of vtkTubeFilters representing the tracts - :type out_list: list - :return: A collection of tracts as a vtkMultiBlockDataSet + :return: The collection of tracts (streamlines) as a vtkMultiBlockDataSet :rtype: vtkMultiBlockDataSet """ + # create a branch and add the streamlines branch = vtk.vtkMultiBlockDataSet() + # create tracts only when at least one was computed # print("Len outlist in root: ", len(out_list)) + # TODO: check if this if statement is required, because we should + # call this function only when tracts exist if not out_list.count(None) == len(out_list): for n, tube in enumerate(out_list): - branch.SetBlock(n, tube.GetOutput()) + branch.SetBlock(n_block + n, tube.GetOutput()) return branch -def tracts_computation(trk_list, root, n_tracts): +def compute_tracts(trk_list, n_tract=0, alpha=255): """Convert the list of all computed tracts given by Trekker run and returns a vtkMultiBlockDataSet :param trk_list: List of lists containing the computed tracts and corresponding coordinates :type trk_list: list - :param root: A collection of tracts as a vtkMultiBlockDataSet - :type root: vtkMultiBlockDataSet - :param n_tracts: - :type n_tracts: int + :param n_tract: The integer ID of the block in the vtkMultiBlockDataSet + :type n_tract: int + :param alpha: The transparency of the streamlines from 0 to 255 (transparent to opaque) + :type alpha: int :return: The updated collection of tracts as a vtkMultiBlockDataSet :rtype: vtkMultiBlockDataSet """ # Transform tracts to array trk_arr = [np.asarray(trk_n).T if trk_n else None for trk_n in trk_list] - # Compute the directions - trk_dir = [compute_directions(trk_n) for trk_n in trk_arr] - + trk_dir = [compute_directions(trk_n, alpha) for trk_n in trk_arr] # Compute the vtk tubes out_list = [compute_tubes(trk_arr_n, trk_dir_n) for trk_arr_n, trk_dir_n in zip(trk_arr, trk_dir)] + # create a branch and add the tracts + branch = create_branch(out_list, n_tract) - root = combine_tracts_root(out_list, root, n_tracts) - - return root + return branch -def compute_tracts(trekker, position, affine, affine_vtk, n_tracts): +def compute_and_visualize_tracts(trekker, position, affine, affine_vtk, n_tracts_max): """ Compute tractograms using the Trekker library. :param trekker: Trekker library instance @@ -183,52 +162,51 @@ def compute_tracts(trekker, position, affine, affine_vtk, n_tracts): :type affine: numpy.ndarray :param affine_vtk: vtkMatrix4x4 isntance with affine transformation matrix :type affine_vtk: vtkMatrix4x4 - :param n_tracts: number of tracts to compute - :type n_tracts: int + :param n_tracts_max: maximum number of tracts to compute + :type n_tracts_max: int """ - # during neuronavigation, root needs to be initialized outside the while loop so the new tracts - # can be appended to the root block set - root = vtk.vtkMultiBlockDataSet() + # root = vtk.vtkMultiBlockDataSet() # Juuso's # seed = np.array([[-8.49, -8.39, 2.5]]) # Baran M1 # seed = np.array([[27.53, -77.37, 46.42]]) seed_trk = img_utils.convert_world_to_voxel(position, affine) - # print("seed example: {}".format(seed_trk)) - trekker.seed_coordinates(np.repeat(seed_trk, n_tracts, axis=0)) - # print("trk list len: ", len(trekker.run())) - trk_list = trekker.run() - if trk_list: - root = tracts_computation(trk_list, root, 0) - Publisher.sendMessage('Remove tracts') - Publisher.sendMessage('Update tracts', root=root, affine_vtk=affine_vtk, coord_offset=position) - else: - Publisher.sendMessage('Remove tracts') - - -def tracts_computation_branch(trk_list, alpha=255): - """Convert the list of all computed tracts given by Trekker run and returns a vtkMultiBlockDataSet - - :param trk_list: List of lists containing the computed tracts and corresponding coordinates - :type trk_list: list - :param alpha: opacity value in the interval [0, 255]. The 0 is no opacity (total transparency). - :type alpha: int - :return: The collection of tracts as a vtkMultiBlockDataSet - :rtype: vtkMultiBlockDataSet - """ - # Transform tracts to array - trk_arr = [np.asarray(trk_n).T if trk_n else None for trk_n in trk_list] - # Compute the directions - trk_dir = [compute_directions(trk_n, alpha) for trk_n in trk_arr] - # Compute the vtk tubes - tube_list = [compute_tubes(trk_arr_n, trk_dir_n) for trk_arr_n, trk_dir_n in zip(trk_arr, trk_dir)] - branch = combine_tracts_branch(tube_list) - - return branch + bundle = vtk.vtkMultiBlockDataSet() + n_branches, n_tracts, count_loop = 0, 0, 0 + n_threads = 2 * const.N_CPU - 1 + + while n_tracts < n_tracts_max: + n_param = 1 + (count_loop % 10) + # rescale the alpha value that defines the opacity of the branch + # the n interval is [1, 10] and the new interval is [51, 255] + # the new interval is defined to have no 0 opacity (minimum is 51, i.e., 20%) + alpha = (n_param - 1) * (255 - 51) / (10 - 1) + 51 + trekker.minFODamp(n_param * 0.01) + + # print("seed example: {}".format(seed_trk)) + trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0)) + # print("trk list len: ", len(trekker.run())) + trk_list = trekker.run() + n_tracts += len(trk_list) + if len(trk_list): + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha) + bundle.SetBlock(n_branches, branch) + n_branches += 1 + + count_loop += 1 + + if (count_loop == 20) and (n_tracts == 0): + break + + Publisher.sendMessage('Remove tracts') + if n_tracts: + Publisher.sendMessage('Update tracts', root=bundle, affine_vtk=affine_vtk, + coord_offset=position, coord_offset_w=seed_trk[0].tolist()) class ComputeTractsThread(threading.Thread): + # TODO: Remove this class and create a case where no ACT is provided in the class ComputeTractsACTThread def __init__(self, inp, queues, event, sle): """Class (threading) to compute real time tractography data for visualization. @@ -265,6 +243,7 @@ class ComputeTractsThread(threading.Thread): trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.inp # n_threads = n_tracts_total + n_threads = int(n_threads/4) p_old = np.array([[0., 0., 0.]]) n_tracts = 0 @@ -312,13 +291,13 @@ class ComputeTractsThread(threading.Thread): # run the trekker, this is the slowest line of code, be careful to just use once! trk_list = trekker.run() - if trk_list: + if len(trk_list) > 2: # print("dist: {}".format(dist)) if dist >= seed_radius: # when moving the coil further than the seed_radius restart the bundle computation bundle = vtk.vtkMultiBlockDataSet() n_branches = 0 - branch = tracts_computation_branch(trk_list) + branch = compute_tracts(trk_list, n_tract=0, alpha=255) bundle.SetBlock(n_branches, branch) n_branches += 1 n_tracts = branch.GetNumberOfBlocks() @@ -326,7 +305,7 @@ class ComputeTractsThread(threading.Thread): # TODO: maybe keep computing even if reaches the maximum elif dist < seed_radius and n_tracts < n_tracts_total: # compute tracts blocks and add to bungle until reaches the maximum number of tracts - branch = tracts_computation_branch(trk_list) + branch = compute_tracts(trk_list, n_tract=0, alpha=255) if bundle: bundle.SetBlock(n_branches, branch) n_tracts += branch.GetNumberOfBlocks() @@ -361,314 +340,218 @@ class ComputeTractsThread(threading.Thread): class ComputeTractsACTThread(threading.Thread): - def __init__(self, inp, queues, event, sle): + def __init__(self, input_list, queues, event, sleep_thread): """Class (threading) to compute real time tractography data for visualization. Tracts are computed using the Trekker library by Baran Aydogan (https://dmritrekker.github.io/) For VTK visualization, each tract (fiber) is a constructed as a tube and many tubes combined in one vtkMultiBlockDataSet named as a branch. Several branches are combined in another vtkMultiBlockDataSet named as bundle, to obtain fast computation and visualization. The bundle dataset is mapped to a single vtkActor. - Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the invesalius 3D scene. + Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the + invesalius 3D scene. Sleep function in run method is used to avoid blocking GUI and more fluent, real-time navigation - :param inp: List of inputs: trekker instance, affine numpy array, seed_offset, seed_radius, n_threads - :type inp: list + :param input_list: List of inputs: trekker instance, affine numpy array, seed offset, total number of tracts, + seed radius, number of threads in computer, ACT data array, affine vtk matrix, + image shift for vtk to mri transformation + :type input_list: list :param queues: Queue list with coord_tracts_queue (Queue instance that manage co-registered coordinates) and tracts_queue (Queue instance that manage the tracts to be visualized) :type queues: list[queue.Queue, queue.Queue] :param event: Threading event to coordinate when tasks as done and allow UI release :type event: threading.Event - :param sle: Sleep pause in seconds - :type sle: float + :param sleep_thread: Sleep pause in seconds + :type sleep_thread: float """ threading.Thread.__init__(self, name='ComputeTractsThreadACT') - self.inp = inp - # self.coord_queue = coord_queue + self.input_list = input_list self.coord_tracts_queue = queues[0] self.tracts_queue = queues[1] - # on first pilots (january 12, 2021) used (-4, 4) - self.coord_list_w = img_utils.create_grid((-2, 2), (0, 20), inp[2]-5, 1) - # self.coord_list_sph = img_utils.create_spherical_grid(10, 1) - # self.coord_list_sph = img_utils.create_spherical_grid(10, 1) - # x_norm = np.linspace(norm.ppf(0.01), norm.ppf(0.99), 2*self.coord_list_sph.shape[0]) - # self.pdf = np.flipud(norm.pdf(x_norm[:self.coord_list_sph.shape[0]], loc=0, scale=2.)) - # self.sph_idx = np.linspace(0, self.coord_list_sph.shape[0], num=self.coord_list_sph.shape[0], dtype=int) - # self.visualization_queue = visualization_queue self.event = event - self.sle = sle - - # prj_data = prj.Project() - # matrix_shape = tuple(prj_data.matrix_shape) - # self.img_shift = matrix_shape[1] + self.sleep_thread = sleep_thread def run(self): - # trekker, affine, offset, n_tracts_total, seed_radius, n_threads = self.inp - trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.inp + trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.input_list - # n_threads = n_tracts_total p_old = np.array([[0., 0., 0.]]) - p_old_pre = np.array([[0., 0., 0.]]) - coord_offset = None - n_tracts = 0 - n_branches = 0 + n_branches, n_tracts, count_loop = 0, 0, 0 bundle = None - sph_sampling = True dist_radius = 1.5 - xyz = img_utils.random_sample_sphere(radius=seed_radius, size=100) - coord_list_sph = np.hstack([xyz, np.ones([xyz.shape[0], 1])]).T + # TODO: Try a denser and bigger grid, because it's just a matrix multiplication + # maybe 15 mm below the coil offset by default and 5 cm deep + # create the rectangular grid to find the gray-white matter boundary + coord_list_w = img_utils.create_grid((-2, 2), (0, 20), offset - 5, 1) + + # create the spherical grid to sample around the seed location + samples_in_sphere = img_utils.random_sample_sphere(radius=seed_radius, size=100) + coord_list_sphere = np.hstack([samples_in_sphere, np.ones([samples_in_sphere.shape[0], 1])]).T m_seed = np.identity(4) + # Compute the tracts - # print('ComputeTractsThread: event {}'.format(self.event.is_set())) while not self.event.is_set(): try: - # print("Computing tracts") # get from the queue the coordinates, coregistration transformation matrix, and flipped matrix m_img_flip = self.coord_tracts_queue.get_nowait() - # coord, m_img, m_img_flip = self.coord_queue.get_nowait() - # print('ComputeTractsThread: get {}'.format(count)) - - # TODO: Remove this is not needed - # 20200402: in this new refactored version the m_img comes different than the position - # the new version m_img is already flixped in y, which means that Y is negative - # if only the Y is negative maybe no need for the flip_x funtcion at all in the navigation - # but check all coord_queue before why now the m_img comes different than position - # 20200403: indeed flip_x is just a -1 multiplication to the Y coordinate, remove function flip_x - # m_img_flip = m_img.copy() - # m_img_flip[1, -1] = -m_img_flip[1, -1] # DEBUG: Uncomment the m_img_flip below so that distance is fixed and tracts keep computing # m_img_flip[:3, -1] = (5., 10., 12.) - dist = abs(np.linalg.norm(p_old_pre - np.asarray(m_img_flip[:3, -1]))) - p_old_pre = m_img_flip[:3, -1].copy() - - # Uncertanity visualization -- - # each tract branch is computed with one set of parameters ajusted from 1 to 10 - n_param = 1 + (n_branches % 10) + dist = abs(np.linalg.norm(p_old - np.asarray(m_img_flip[:3, -1]))) + p_old = m_img_flip[:3, -1].copy() + + # Uncertainty visualization -- + # each tract branch is computed with one minFODamp adjusted from 0.01 to 0.1 + # the higher the minFODamp the more the streamlines are faithful to the data, so they become more strict + # but also may loose some connections. + # the lower the more relaxed streamline also increases the chance of false positives + n_param = 1 + (count_loop % 10) # rescale the alpha value that defines the opacity of the branch # the n interval is [1, 10] and the new interval is [51, 255] # the new interval is defined to have no 0 opacity (minimum is 51, i.e., 20%) alpha = (n_param - 1) * (255 - 51) / (10 - 1) + 51 trekker.minFODamp(n_param * 0.01) - trekker.dataSupportExponent(n_param * 0.1) # --- + try: + # The original seed location is replaced by the gray-white matter interface that is closest to + # the coil center + coord_list_w_tr = m_img_flip @ coord_list_w + coord_offset = grid_offset(act_data, coord_list_w_tr, img_shift) + except IndexError: + # This error might be caused by the coordinate exceeding the image array dimensions. + # This happens during navigation because the coil location can have coordinates outside the image + # boundaries + # Translate the coordinate along the normal vector of the object/coil + coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2] + # --- + + # Spherical sampling of seed coordinates --- + # compute the samples of a sphere centered on seed coordinate offset by the grid + # given in the invesalius-vtk space + samples = np.random.choice(coord_list_sphere.shape[1], size=100) + m_seed[:-1, -1] = coord_offset.copy() + # translate the spherical grid samples to the coil location in invesalius-vtk space + seed_trk_r_inv = m_seed @ coord_list_sphere[:, samples] + + coord_offset_w = np.hstack((coord_offset, 1.0)).reshape([4, 1]) + + try: + # Anatomic constrained seed computation --- + # find only the samples inside the white matter as with the ACT enabled in the Trekker, + # it does not compute any tracts outside the white matter + # convert from inveslaius-vtk to mri space + seed_trk_r_mri = seed_trk_r_inv[:3, :].T.astype(int) + np.array([[0, img_shift, 0]], dtype=np.int32) + labs = act_data[seed_trk_r_mri[..., 0], seed_trk_r_mri[..., 1], seed_trk_r_mri[..., 2]] + # find all samples in the white matter + labs_id = np.where(labs == 1) + # Trekker has internal multiprocessing approach done in C. Here the number of available threads - 1 + # is given, but in case a large number of tracts is requested, it will compute all in parallel + # automatically for a more fluent navigation, better to compute the maximum number the computer can + # handle otherwise gets too slow for the multithreading in Python + seed_trk_r_inv_sampled = seed_trk_r_inv[:, labs_id[0][:n_threads]] + + except IndexError: + # same as on the grid offset above, if the coil is too far from the mri volume the array indices + # are beyond the mri boundaries + # in this case use the grid center instead of the spherical samples + seed_trk_r_inv_sampled = coord_offset_w.copy() + + # convert to the world coordinate system for trekker + seed_trk_r_world_sampled = np.linalg.inv(affine) @ seed_trk_r_inv_sampled + seed_trk_r_world_sampled = seed_trk_r_world_sampled.T[:, :3] + + # convert to the world coordinate system for saving in the marker list + coord_offset_w = np.linalg.inv(affine) @ coord_offset_w + coord_offset_w = np.squeeze(coord_offset_w.T[:, :3]) + + # DEBUG: uncomment the seed_trk below + # seed_trk.shape == [1, 3] + # Juuso's + # seed_trk = np.array([[-8.49, -8.39, 2.5]]) + # Baran M1 + # seed_trk = np.array([[27.53, -77.37, 46.42]]) + # print("Given: {}".format(seed_trk.shape)) + # print("Seed: {}".format(seed)) + # joonas seed that has good tracts + # seed_trk = np.array([[29.12, -13.33, 31.65]]) + # seed_trk_img = np.array([[117, 127, 161]]) + # When moving the coil further than the seed_radius restart the bundle computation # Currently, it stops to compute tracts when the maximum number of tracts is reached maybe keep # computing even if reaches the maximum if dist >= dist_radius: - # Anatomic constrained seed computation --- - # The original seed location is replaced by the gray-white matter interface that is closest to - # the coil center - try: - #TODO: Create a dialog error to say when the ACT data is not loaded and prevent - # the interface from freezing. Give the user a chance to load it (maybe in task_navigator) - coord_list_w_tr = m_img_flip @ self.coord_list_w - coord_offset = grid_offset(act_data, coord_list_w_tr, img_shift) - except IndexError: - # This error might be caused by the coordinate exceeding the image array dimensions. - # Needs further verification. - coord_offset = None - # --- - - # Translate the coordinate along the normal vector of the object/coil --- - if coord_offset is None: - # apply the coil transformation matrix - coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2] - # --- - - # convert the world coordinates to the voxel space for using as a seed in Trekker - # seed_trk.shape == [1, 3] - seed_trk = img_utils.convert_world_to_voxel(coord_offset, affine) - # print("Desired: {}".format(seed_trk.shape)) - - # DEBUG: uncomment the seed_trk below - # Juuso's - # seed_trk = np.array([[-8.49, -8.39, 2.5]]) - # Baran M1 - # seed_trk = np.array([[27.53, -77.37, 46.42]]) - # print("Given: {}".format(seed_trk.shape)) - # print("Seed: {}".format(seed)) - - # Spherical sampling of seed coordinates --- - if sph_sampling: - # CHECK: We use ACT only for the origin seed, but not for all the other coordinates. - # Check how this can be solved. Applying ACT to all coordinates is maybe too much. - # Maybe it doesn't matter because the ACT is just to help finding the closest location to - # the TMS coil center. Also, note that the spherical sampling is applied only when the coil - # location changes, all further iterations used the fixed n_threads samples to compute the - # remaining tracts. - - # samples = np.random.choice(self.sph_idx, size=n_threads, p=self.pdf) - # m_seed[:-1, -1] = seed_trk - # sph_seed = m_seed @ self.coord_list_sph - # seed_trk_r = sph_seed[samples, :] - samples = np.random.choice(coord_list_sph.shape[1], size=n_threads) - m_seed[:-1, -1] = seed_trk - seed_trk_r = m_seed @ coord_list_sph[:, samples] - seed_trk_r = seed_trk_r[:-1, :].T - else: - # set the seeds for trekker, one seed is repeated n_threads times - # shape (24, 3) - seed_trk_r = np.repeat(seed_trk, n_threads, axis=0) - - # --- - - # Trekker has internal multiprocessing approach done in C. Here the number of available threads is - # given, but in case a large number of tracts is requested, it will compute all in parallel - # automatically for a more fluent navigation, better to compute the maximum number the computer - # handles - trekker.seed_coordinates(seed_trk_r) - # trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0)) + bundle = None + n_tracts, n_branches = 0, 0 + # we noticed that usually the navigation lags or crashes when moving the coil location + # to reduce the overhead for when the coil is moving, we compute only half the number of tracts + # that should be computed when the coil is fixed in the same location + # required input is Nx3 array + trekker.seed_coordinates(seed_trk_r_world_sampled[::2, :]) # run the trekker, this is the slowest line of code, be careful to just use once! trk_list = trekker.run() # check if any tract was found, otherwise doesn't count - if trk_list: + if len(trk_list): + # a bundle consists for multiple branches and each branch consists of multiple streamlines + # every iteration in the main loop adds a branch to the bundle + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha) + n_tracts = branch.GetNumberOfBlocks() + + # create and add branch to the bundle bundle = vtk.vtkMultiBlockDataSet() - branch = tracts_computation_branch(trk_list, alpha) bundle.SetBlock(n_branches, branch) n_branches = 1 - n_tracts = branch.GetNumberOfBlocks() - else: - bundle = None - n_branches = 0 - n_tracts = 0 elif dist < dist_radius and n_tracts < n_tracts_total: + # when the coil is fixed in place and the number of tracts is smaller than the total + if not bundle: + # same as above, when creating the bundle (vtkMultiBlockDataSet) we only compute half the number + # of tracts to reduce the overhead + bundle = vtk.vtkMultiBlockDataSet() + # required input is Nx3 array + trekker.seed_coordinates(seed_trk_r_world_sampled[::2, :]) + n_tracts, n_branches = 0, 0 + else: + # if the bundle exists compute all tracts requested + # required input is Nx3 array + trekker.seed_coordinates(seed_trk_r_world_sampled) + trk_list = trekker.run() - if trk_list: + + if len(trk_list): # compute tract blocks and add to bundle until reaches the maximum number of tracts # the alpha changes depending on the parameter set - branch = tracts_computation_branch(trk_list, alpha) - if bundle: - bundle.SetBlock(n_branches, branch) - n_tracts += branch.GetNumberOfBlocks() - n_branches += 1 - else: - bundle = vtk.vtkMultiBlockDataSet() - bundle.SetBlock(n_branches, branch) - n_branches = 1 - n_tracts = branch.GetNumberOfBlocks() - # else: - # bundle = None + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha) + n_tracts += branch.GetNumberOfBlocks() + # add branch to the bundle + bundle.SetBlock(n_branches, branch) + n_branches += 1 - # else: - # bundle = None + # keep adding to the number of loops even if the tracts were not find + # this will keep the minFODamp changing and new seed coordinates being tried which would allow + # higher probability of finding tracts + count_loop += 1 - # rethink if this should be inside the if condition, it may lock the thread if no tracts are found - # use no wait to ensure maximum speed and avoid visualizing old tracts in the queue, this might + # use "nowait" to ensure maximum speed and avoid visualizing old tracts in the queue, this might # be more evident in slow computer or for heavier tract computations, it is better slow update # than visualizing old data - # self.visualization_queue.put_nowait([coord, m_img, bundle]) - self.tracts_queue.put_nowait((bundle, affine_vtk, coord_offset)) - # print('ComputeTractsThread: put {}'.format(count)) - + self.tracts_queue.put_nowait((bundle, affine_vtk, coord_offset, coord_offset_w)) self.coord_tracts_queue.task_done() - # self.coord_queue.task_done() - # print('ComputeTractsThread: done {}'.format(count)) # sleep required to prevent user interface from being unresponsive - time.sleep(self.sle) + time.sleep(self.sleep_thread) # if no coordinates pass except queue.Empty: - # print("Empty queue in tractography") pass # if queue is full mark as done (may not be needed in this new "nowait" method) except queue.Full: - # self.coord_queue.task_done() self.coord_tracts_queue.task_done() -class ComputeTractsThreadSingleBlock(threading.Thread): - - def __init__(self, inp, affine_vtk, coord_queue, visualization_queue, event, sle): - """Class (threading) to compute real time tractography data for visualization in a single loop. - - Different than ComputeTractsThread because it does not keep adding tracts to the bundle until maximum, - is reached. It actually compute all requested tracts at once. (Might be deleted in the future)! - Tracts are computed using the Trekker library by Baran Aydogan (https://dmritrekker.github.io/) - For VTK visualization, each tract (fiber) is a constructed as a tube and many tubes combined in one - vtkMultiBlockDataSet named as a branch. Several branches are combined in another vtkMultiBlockDataSet named as - bundle, to obtain fast computation and visualization. The bundle dataset is mapped to a single vtkActor. - Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the invesalius 3D scene. - - Sleep function in run method is used to avoid blocking GUI and more fluent, real-time navigation - - :param inp: List of inputs: trekker instance, affine numpy array, seed_offset, seed_radius, n_threads - :type inp: list - :param affine_vtk: Affine matrix in vtkMatrix4x4 instance to update objects position in 3D scene - :type affine_vtk: vtkMatrix4x4 - :param coord_queue: Queue instance that manage coordinates read from tracking device and coregistered - :type coord_queue: queue.Queue - :param visualization_queue: Queue instance that manage coordinates to be visualized - :type visualization_queue: queue.Queue - :param event: Threading event to coordinate when tasks as done and allow UI release - :type event: threading.Event - :param sle: Sleep pause in seconds - :type sle: float - """ - - threading.Thread.__init__(self, name='ComputeTractsThread') - self.inp = inp - self.affine_vtk = affine_vtk - self.coord_queue = coord_queue - self.visualization_queue = visualization_queue - self.event = event - self.sle = sle - - def run(self): - - trekker, affine, offset, n_tracts_total, seed_radius, n_threads = self.inp - # as a single block, computes all the maximum number of tracts at once, not optimal for navigation - n_threads = n_tracts_total - p_old = np.array([[0., 0., 0.]]) - root = vtk.vtkMultiBlockDataSet() - - # Compute the tracts - # print('ComputeTractsThread: event {}'.format(self.event.is_set())) - while not self.event.is_set(): - try: - coord, [coord_raw, markers_flag], m_img, m_img_flip = self.coord_queue.get_nowait() - - # translate the coordinate along the normal vector of the object/coil - coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2] - # coord_offset = np.array([[27.53, -77.37, 46.42]]) - dist = abs(np.linalg.norm(p_old - np.asarray(coord_offset))) - p_old = coord_offset.copy() - seed_trk = img_utils.convert_world_to_voxel(coord_offset, affine) - # Juuso's - # seed_trk = np.array([[-8.49, -8.39, 2.5]]) - # Baran M1 - # seed_trk = np.array([[27.53, -77.37, 46.42]]) - # print("Seed: {}".format(seed)) - - # set the seeds for trekker, one seed is repeated n_threads times - # trekker has internal multiprocessing approach done in C. Here the number of available threads is give, - # but in case a large number of tracts is requested, it will compute all in parallel automatically - # for a more fluent navigation, better to compute the maximum number the computer handles - trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0)) - # run the trekker, this is the slowest line of code, be careful to just use once! - trk_list = trekker.run() - - if trk_list: - # if the seed is outside the defined radius, restart the bundle computation - if dist >= seed_radius: - root = tracts_computation(trk_list, root, 0) - self.visualization_queue.put_nowait((coord, m_img, root)) - - self.coord_queue.task_done() - time.sleep(self.sle) - except queue.Empty: - pass - except queue.Full: - self.coord_queue.task_done() - - def set_trekker_parameters(trekker, params): """Set all user-defined parameters for tractography computation using the Trekker library @@ -680,21 +563,27 @@ def set_trekker_parameters(trekker, params): :rtype: list """ trekker.seed_maxTrials(params['seed_max']) - # trekker.stepSize(params['step_size']) - trekker.minFODamp(params['min_fod']) - # trekker.probeQuality(params['probe_quality']) - # trekker.maxEstInterval(params['max_interval']) - # trekker.minRadiusOfCurvature(params['min_radius_curv']) - # trekker.probeLength(params['probe_length']) + trekker.stepSize(params['step_size']) + # minFODamp is not set because it should vary in the loop to create the + # different transparency tracts + # trekker.minFODamp(params['min_fod']) + trekker.probeQuality(params['probe_quality']) + trekker.maxEstInterval(params['max_interval']) + trekker.minRadiusOfCurvature(params['min_radius_curvature']) + trekker.probeLength(params['probe_length']) trekker.writeInterval(params['write_interval']) - trekker.maxLength(params['max_lenth']) - trekker.minLength(params['min_lenth']) + # these two does not need to be set in the new package + # trekker.maxLength(params['max_length']) + trekker.minLength(params['min_length']) trekker.maxSamplingPerStep(params['max_sampling_step']) + trekker.dataSupportExponent(params['data_support_exponent']) + #trekker.useBestAtInit(params['use_best_init']) + #trekker.initMaxEstTrials(params['init_max_est_trials']) # check number if number of cores is valid in configuration file, # otherwise use the maximum number of threads which is usually 2*N_CPUS - n_threads = 2 * const.N_CPU - if isinstance((params['numb_threads']), int) and params['numb_threads'] <= 2*const.N_CPU: + n_threads = 2 * const.N_CPU - 1 + if isinstance((params['numb_threads']), int) and params['numb_threads'] <= (2*const.N_CPU-1): n_threads = params['numb_threads'] trekker.numberOfThreads(n_threads) @@ -712,13 +601,21 @@ def grid_offset(data, coord_list_w_tr, img_shift): # extract the first occurrence of a specific label from the MRI image labs = data[coord_list_w_tr_mri[..., 0], coord_list_w_tr_mri[..., 1], coord_list_w_tr_mri[..., 2]] - lab_first = np.argmax(labs == 1) - if labs[lab_first] == 1: - pt_found = coord_list_w_tr_mri[lab_first, :] + lab_first = np.where(labs == 1) + if not lab_first: + pt_found_inv = None + else: + pt_found = coord_list_w_tr[:, lab_first[0][0]][:3] # convert coordinate back to invesalius 3D space pt_found_inv = pt_found - np.array([0., img_shift, 0.]) - else: - pt_found_inv = None + + # lab_first = np.argmax(labs == 1) + # if labs[lab_first] == 1: + # pt_found = coord_list_w_tr_mri[lab_first, :] + # # convert coordinate back to invesalius 3D space + # pt_found_inv = pt_found - np.array([0., img_shift, 0.]) + # else: + # pt_found_inv = None # # convert to world coordinate space to use as seed for fiber tracking # pt_found_tr = np.append(pt_found, 1)[np.newaxis, :].T diff --git a/invesalius/data/viewer_volume.py b/invesalius/data/viewer_volume.py index 0e97b4a..aa75666 100644 --- a/invesalius/data/viewer_volume.py +++ b/invesalius/data/viewer_volume.py @@ -307,7 +307,6 @@ class Viewer(wx.Panel): Publisher.subscribe(self.OnRemoveTracts, 'Remove tracts') Publisher.subscribe(self.UpdateSeedOffset, 'Update seed offset') Publisher.subscribe(self.UpdateMarkerOffsetState, 'Update marker offset state') - Publisher.subscribe(self.UpdateMarkerOffsetPosition, 'Update marker offset') Publisher.subscribe(self.AddPeeledSurface, 'Update peel') Publisher.subscribe(self.GetPeelCenters, 'Get peel centers and normals') Publisher.subscribe(self.Initlocator_viewer, 'Get init locator') @@ -857,8 +856,10 @@ class Viewer(wx.Panel): else: self.DisableCoilTracker() if self.actor_peel: - self.object_orientation_torus_actor.SetVisibility(1) - self.obj_projection_arrow_actor.SetVisibility(1) + if self.object_orientation_torus_actor: + self.object_orientation_torus_actor.SetVisibility(1) + if self.obj_projection_arrow_actor: + self.obj_projection_arrow_actor.SetVisibility(1) def OnUpdateObjectTargetGuide(self, m_img, coord): @@ -1607,10 +1608,6 @@ class Viewer(wx.Panel): self.ren.AddActor(self.mark_actor) self.Refresh() - def UpdateMarkerOffsetPosition(self, coord_offset): - self.mark_actor.SetPosition(coord_offset) - self.Refresh() - def UpdateObjectOrientation(self, m_img, coord): # print("Update object orientation") @@ -1687,7 +1684,7 @@ class Viewer(wx.Panel): # self.ball_actor.SetVisibility(1) self.Refresh() - def OnUpdateTracts(self, root=None, affine_vtk=None, coord_offset=None): + def OnUpdateTracts(self, root=None, affine_vtk=None, coord_offset=None, coord_offset_w=None): mapper = vtk.vtkCompositePolyDataMapper2() mapper.SetInputDataObject(root) diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index db7c8cb..5f2471e 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -4203,7 +4203,7 @@ class GoToDialogScannerCoord(wx.Dialog): main_sizer.Fit(self) self.orientation = None - self.affine = None + self.affine = np.identity(4) self.__bind_events() diff --git a/invesalius/gui/task_navigator.py b/invesalius/gui/task_navigator.py index e063292..6b9408a 100644 --- a/invesalius/gui/task_navigator.py +++ b/invesalius/gui/task_navigator.py @@ -165,7 +165,7 @@ class InnerFoldPanel(wx.Panel): # Study this. fold_panel = fpb.FoldPanelBar(self, -1, wx.DefaultPosition, - (10, 310), 0, fpb.FPB_SINGLE_FOLD) + (10, 330), 0, fpb.FPB_SINGLE_FOLD) # Initialize Tracker and PedalConnection objects here so that they are available to several panels. # @@ -1078,17 +1078,15 @@ class ObjectRegistrationPanel(wx.Panel): try: if filename: - #TODO: Improve method to read the file, using "with" similar to OnLoadParameters - data = np.loadtxt(filename, delimiter='\t') - self.obj_fiducials = data[:, :3] - self.obj_orients = data[:, 3:] + with open(filename, 'r') as text_file: + data = [s.split('\t') for s in text_file.readlines()] - text_file = open(filename, "r") - header = text_file.readline().split('\t') - text_file.close() + registration_coordinates = np.array(data[1:]).astype(np.float32) + self.obj_fiducials = registration_coordinates[:, :3] + self.obj_orients = registration_coordinates[:, 3:] - self.obj_name = header[1] - self.obj_ref_mode = int(header[-1]) + self.obj_name = data[0][1] + self.obj_ref_mode = int(data[0][-1]) self.checktrack.Enable(1) self.checktrack.SetValue(True) @@ -1098,7 +1096,7 @@ class ObjectRegistrationPanel(wx.Panel): label=_("Object file successfully loaded")) Publisher.sendMessage('Update track object state', flag=True, obj_name=self.obj_name) Publisher.sendMessage('Change camera checkbox', status=False) - # wx.MessageBox(_("Object file successfully loaded"), _("Load")) + wx.MessageBox(_("Object file successfully loaded"), _("InVesalius 3")) except: wx.MessageBox(_("Object registration file incompatible."), _("InVesalius 3")) Publisher.sendMessage('Update status text in GUI', label="") @@ -1449,8 +1447,8 @@ class MarkersPanel(wx.Panel): else: self.nav_status = True - def UpdateSeedCoordinates(self, root=None, affine_vtk=None, coord_offset=(0, 0, 0)): - self.current_seed = coord_offset + def UpdateSeedCoordinates(self, root=None, affine_vtk=None, coord_offset=(0, 0, 0), coord_offset_w=(0, 0, 0)): + self.current_seed = coord_offset_w def UpdateRobotCoordinates(self, coordinates_raw, markers_flag): self.raw_target_robot = coordinates_raw[1], coordinates_raw[2] @@ -1735,7 +1733,7 @@ class TractographyPanel(wx.Panel): default_colour = wx.SystemSettings_GetColour(wx.SYS_COLOUR_MENUBAR) self.SetBackgroundColour(default_colour) - self.affine = None + self.affine = np.identity(4) self.affine_vtk = None self.trekker = None self.n_tracts = const.N_TRACTS @@ -1992,7 +1990,6 @@ class TractographyPanel(wx.Panel): self.nav_status = nav_status def OnLinkBrain(self, event=None): - Publisher.sendMessage('Update status text in GUI', label=_("Busy")) Publisher.sendMessage('Begin busy cursor') mask_path = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, _("Import brain mask")) img_path = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, _("Import T1 anatomical image")) @@ -2006,31 +2003,37 @@ class TractographyPanel(wx.Panel): slic = sl.Slice() prj_data = prj.Project() matrix_shape = tuple(prj_data.matrix_shape) + spacing = tuple(prj_data.spacing) + img_shift = spacing[1] * (matrix_shape[1] - 1) self.affine = slic.affine.copy() - self.affine[1, -1] -= matrix_shape[1] + self.affine[1, -1] -= img_shift self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine) - try: - self.brain_peel = brain.Brain(img_path, mask_path, self.n_peels, self.affine_vtk) - self.brain_actor = self.brain_peel.get_actor(self.peel_depth, self.affine_vtk) - self.brain_actor.GetProperty().SetOpacity(self.brain_opacity) - Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor) - Publisher.sendMessage('Get peel centers and normals', centers=self.brain_peel.peel_centers, - normals=self.brain_peel.peel_normals) - Publisher.sendMessage('Get init locator', locator=self.brain_peel.locator) - self.checkpeeling.Enable(1) - self.checkpeeling.SetValue(True) - self.spin_opacity.Enable(1) - Publisher.sendMessage('Update status text in GUI', label=_("Brain model loaded")) - self.peel_loaded = True - Publisher.sendMessage('Update peel visualization', data= self.peel_loaded) - except: - wx.MessageBox(_("Unable to load brain mask."), _("InVesalius 3")) + if mask_path and img_path: + Publisher.sendMessage('Update status text in GUI', label=_("Busy")) + try: + self.brain_peel = brain.Brain(img_path, mask_path, self.n_peels, self.affine_vtk) + self.brain_actor = self.brain_peel.get_actor(self.peel_depth, self.affine_vtk) + self.brain_actor.GetProperty().SetOpacity(self.brain_opacity) + + self.checkpeeling.Enable(1) + self.checkpeeling.SetValue(True) + self.spin_opacity.Enable(1) + self.peel_loaded = True + + Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor) + Publisher.sendMessage('Get peel centers and normals', centers=self.brain_peel.peel_centers, + normals=self.brain_peel.peel_normals) + Publisher.sendMessage('Get init locator', locator=self.brain_peel.locator) + Publisher.sendMessage('Update status text in GUI', label=_("Brain model loaded")) + Publisher.sendMessage('Update peel visualization', data= self.peel_loaded) + except: + Publisher.sendMessage('Update status text in GUI', label=_("Brain mask initialization failed.")) + wx.MessageBox(_("Unable to load brain mask."), _("InVesalius 3")) Publisher.sendMessage('End busy cursor') def OnLinkFOD(self, event=None): - Publisher.sendMessage('Update status text in GUI', label=_("Busy")) Publisher.sendMessage('Begin busy cursor') filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import Trekker FOD")) # Juuso @@ -2041,69 +2044,83 @@ class TractographyPanel(wx.Panel): # FOD_path = 'Baran_FOD.nii' # filename = os.path.join(data_dir, FOD_path) - # if not self.affine_vtk: - # slic = sl.Slice() - # self.affine = slic.affine - # self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine) - if not self.affine_vtk: slic = sl.Slice() prj_data = prj.Project() matrix_shape = tuple(prj_data.matrix_shape) + spacing = tuple(prj_data.spacing) + img_shift = spacing[1] * (matrix_shape[1] - 1) self.affine = slic.affine.copy() - self.affine[1, -1] -= matrix_shape[1] + self.affine[1, -1] -= img_shift self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine) - # try: - - self.trekker = Trekker.initialize(filename.encode('utf-8')) - self.trekker, n_threads = dti.set_trekker_parameters(self.trekker, self.trekker_cfg) - - self.checktracts.Enable(1) - self.checktracts.SetValue(True) - self.view_tracts = True - Publisher.sendMessage('Update Trekker object', data=self.trekker) - Publisher.sendMessage('Update number of threads', data=n_threads) - Publisher.sendMessage('Update tracts visualization', data=1) - Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialized")) - # except: - # wx.MessageBox(_("Unable to initialize Trekker, check FOD and config files."), _("InVesalius 3")) + if filename: + Publisher.sendMessage('Update status text in GUI', label=_("Busy")) + try: + self.trekker = Trekker.initialize(filename.encode('utf-8')) + self.trekker, n_threads = dti.set_trekker_parameters(self.trekker, self.trekker_cfg) + + self.checktracts.Enable(1) + self.checktracts.SetValue(True) + self.view_tracts = True + + Publisher.sendMessage('Update Trekker object', data=self.trekker) + Publisher.sendMessage('Update number of threads', data=n_threads) + Publisher.sendMessage('Update tracts visualization', data=1) + Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialized")) + # except: + # wx.MessageBox(_("Unable to initialize Trekker, check FOD and config files."), _("InVesalius 3")) + except: + Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialization failed.")) + wx.MessageBox(_("Unable to load FOD."), _("InVesalius 3")) Publisher.sendMessage('End busy cursor') def OnLoadACT(self, event=None): - Publisher.sendMessage('Update status text in GUI', label=_("Busy")) - Publisher.sendMessage('Begin busy cursor') - filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import anatomical labels")) - # Baran - # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\anat_reg_improve_20200609' - # act_path = 'Baran_trekkerACTlabels_inFODspace.nii' - # filename = os.path.join(data_dir, act_path) - - act_data = nb.squeeze_image(nb.load(filename)) - act_data = nb.as_closest_canonical(act_data) - act_data.update_header() - act_data_arr = act_data.get_fdata() + if self.trekker: + Publisher.sendMessage('Begin busy cursor') + filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import anatomical labels")) + # Baran + # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\anat_reg_improve_20200609' + # act_path = 'Baran_trekkerACTlabels_inFODspace.nii' + # filename = os.path.join(data_dir, act_path) + + if not self.affine_vtk: + slic = sl.Slice() + prj_data = prj.Project() + matrix_shape = tuple(prj_data.matrix_shape) + spacing = tuple(prj_data.spacing) + img_shift = spacing[1] * (matrix_shape[1] - 1) + self.affine = slic.affine.copy() + self.affine[1, -1] -= img_shift + self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine) - if not self.affine_vtk: - slic = sl.Slice() - prj_data = prj.Project() - matrix_shape = tuple(prj_data.matrix_shape) - self.affine = slic.affine.copy() - self.affine[1, -1] -= matrix_shape[1] - self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine) - - self.checkACT.Enable(1) - self.checkACT.SetValue(True) - - Publisher.sendMessage('Update ACT data', data=act_data_arr) - Publisher.sendMessage('Enable ACT', data=True) - # Publisher.sendMessage('Create grid', data=act_data_arr, affine=self.affine) - # Publisher.sendMessage('Update number of threads', data=n_threads) - # Publisher.sendMessage('Update tracts visualization', data=1) - Publisher.sendMessage('Update status text in GUI', label=_("Trekker ACT loaded")) - - Publisher.sendMessage('End busy cursor') + try: + Publisher.sendMessage('Update status text in GUI', label=_("Busy")) + if filename: + act_data = nb.squeeze_image(nb.load(filename)) + act_data = nb.as_closest_canonical(act_data) + act_data.update_header() + act_data_arr = act_data.get_fdata() + + self.checkACT.Enable(1) + self.checkACT.SetValue(True) + + # ACT rules should be as follows: + self.trekker.pathway_stop_at_entry(filename.encode('utf-8'), -1) # outside + self.trekker.pathway_discard_if_ends_inside(filename.encode('utf-8'), 1) # wm + self.trekker.pathway_discard_if_enters(filename.encode('utf-8'), 0) # csf + + Publisher.sendMessage('Update ACT data', data=act_data_arr) + Publisher.sendMessage('Enable ACT', data=True) + Publisher.sendMessage('Update status text in GUI', label=_("Trekker ACT loaded")) + except: + Publisher.sendMessage('Update status text in GUI', label=_("ACT initialization failed.")) + wx.MessageBox(_("Unable to load ACT."), _("InVesalius 3")) + + Publisher.sendMessage('End busy cursor') + else: + wx.MessageBox(_("Load FOD image before the ACT."), _("InVesalius 3")) def OnLoadParameters(self, event=None): import json @@ -2146,10 +2163,13 @@ class TractographyPanel(wx.Panel): # print("Running during navigation") coord_flip = list(position[:3]) coord_flip[1] = -coord_flip[1] - dti.compute_tracts(self.trekker, coord_flip, self.affine, self.affine_vtk, - self.n_tracts) + dti.compute_and_visualize_tracts(self.trekker, coord_flip, self.affine, self.affine_vtk, + self.n_tracts) def OnCloseProject(self): + self.trekker = None + self.trekker_cfg = const.TREKKER_CONFIG + self.checktracts.SetValue(False) self.checktracts.Enable(0) self.checkpeeling.SetValue(False) diff --git a/invesalius/navigation/navigation.py b/invesalius/navigation/navigation.py index a6e7142..bc5c299 100644 --- a/invesalius/navigation/navigation.py +++ b/invesalius/navigation/navigation.py @@ -27,7 +27,6 @@ import numpy as np import invesalius.constants as const import invesalius.project as prj import invesalius.data.bases as db -import invesalius.data.coordinates as dco import invesalius.data.coregistration as dcr import invesalius.data.serial_port_connection as spc import invesalius.data.slice_ as sl @@ -98,12 +97,11 @@ class UpdateNavigationScene(threading.Thread): # use of CallAfter is mandatory otherwise crashes the wx interface if self.view_tracts: - bundle, affine_vtk, coord_offset = self.tracts_queue.get_nowait() + bundle, affine_vtk, coord_offset, coord_offset_w = self.tracts_queue.get_nowait() #TODO: Check if possible to combine the Remove tracts with Update tracts in a single command wx.CallAfter(Publisher.sendMessage, 'Remove tracts') - wx.CallAfter(Publisher.sendMessage, 'Update tracts', root=bundle, - affine_vtk=affine_vtk, coord_offset=coord_offset) - # wx.CallAfter(Publisher.sendMessage, 'Update marker offset', coord_offset=coord_offset) + wx.CallAfter(Publisher.sendMessage, 'Update tracts', root=bundle, affine_vtk=affine_vtk, + coord_offset=coord_offset, coord_offset_w=coord_offset_w) self.tracts_queue.task_done() if self.serial_port_enabled: @@ -118,10 +116,11 @@ class UpdateNavigationScene(threading.Thread): wx.CallAfter(Publisher.sendMessage, 'Set cross focal point', position=coord) wx.CallAfter(Publisher.sendMessage, 'Update raw coordinates', coordinates_raw=coordinates_raw, markers_flag=markers_flag) wx.CallAfter(Publisher.sendMessage, 'Update slice viewer') + wx.CallAfter(Publisher.sendMessage, 'Sensor ID', markers_flag=markers_flag) if view_obj: wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) - wx.CallAfter(Publisher.sendMessage, 'Update object arrow matrix',m_img=m_img, coord=coord, flag= self.peel_loaded) + wx.CallAfter(Publisher.sendMessage, 'Update object arrow matrix', m_img=m_img, coord=coord, flag= self.peel_loaded) self.coord_queue.task_done() # print('UpdateScene: done {}'.format(count)) # count += 1 @@ -191,7 +190,7 @@ class Navigation(): def UpdateSleep(self, sleep): self.sleep_nav = sleep - self.serial_port_connection.sleep_nav = sleep + # self.serial_port_connection.sleep_nav = sleep def UpdateSerialPort(self, serial_port_in_use, com_port=None, baud_rate=None): self.serial_port_in_use = serial_port_in_use @@ -307,15 +306,22 @@ class Navigation(): if self.view_tracts: # initialize Trekker parameters + # TODO: This affine and affine_vtk are created 4 times. To improve, create a new affine object inside + # Slice() that contains the transformation with the img_shift. Rename it to avoid confusion to the + # affine, for instance it can be: affine_world_to_invesalius_vtk slic = sl.Slice() prj_data = prj.Project() matrix_shape = tuple(prj_data.matrix_shape) + spacing = tuple(prj_data.spacing) + img_shift = spacing[1] * (matrix_shape[1] - 1) affine = slic.affine.copy() - affine[1, -1] -= matrix_shape[1] + affine[1, -1] -= img_shift affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(affine) + Publisher.sendMessage("Update marker offset state", create=True) + self.trk_inp = self.trekker, affine, self.seed_offset, self.n_tracts, self.seed_radius,\ - self.n_threads, self.act_data, affine_vtk, matrix_shape[1] + self.n_threads, self.act_data, affine_vtk, img_shift # print("Appending the tract computation thread!") queues = [self.coord_tracts_queue, self.tracts_queue] if self.enable_act: -- libgit2 0.21.2