Download as
Browse Code »

Commit 46e9180122ff30cd6e3d4c95c4a597e2fdd1ffab

Authored by Victor Hugo Souza
Committed by GitHub
1 parent 3b8f8e1f
Exists in master

ENH: Update tractography parameters and fix config file load (#400) 

* ENH: Update tractography parameters and fix config file load

* FIX: Problem when loading the object registration after loading FOD

* FIX: Fix torus bug when creating surface and update tractography parameters

* FIX: Include spacing in the offset for tract computation

* ADD: Update tractography parameters and image shift with spacing

* ENH: Improve grid offset for tractography

* WIP: Improved tractography computation and seed sampling

* WIP: Refactored and cleaned ACT tractography needs testing in navigation

* FIX: Bump conda environment packages to python 3.8 and tract parameters

* FIX: Save the seed marker in world coordinate space

* ENH: Affine is identity by default and not None

Co-authored-by: Renan <renan_hiroshi@hotmail.com>
Showing 12 changed files with 401 additions and 450 deletions   Show diff stats
environment.yml
  Diff comments   View file @46e9180
... ... @@ -2,7 +2,7 @@ channels:
2 2 - conda-forge
3 3 - bioconda
4 4 dependencies:
5   - - python=3.7
  5 + - python=3.8
6 6 - cython==0.29.24
7 7 - pillow==8.3.2
8 8 - pypubsub==4.0.3
... ... @@ -17,13 +17,21 @@ dependencies:
17 17 - scipy==1.7.1
18 18 - vtk==9.0.3
19 19 - wxpython==4.1.1
  20 + - mido==1.2.10
  21 + - nest-asyncio==1.5.1
20 22 - pip
21 23 - pip:
22   - - python-gdcm==3.0.9.1
  24 + - aioconsole==0.3.2
  25 + - opencv-python==4.5.3.56
23 26 - plaidml-keras==0.7.0
24   - - theano==1.0.5
25   - - pyacvd==0.2.3
  27 + - python-gdcm==3.0.9.1
  28 + - python-rtmidi==1.4.9
  29 + - python-socketio[client]==5.3.0
  30 + - pyacvd==0.2.7
26 31 - pyclaron
27 32 - polhemusFT
28 33 - polhemus
29 34 - pypolaris
  35 + - theano==1.0.5
  36 + - uvicorn[standard]==0.15.0
  37 + - https://github.com/baranaydogan/trekker/raw/master/binaries/Trekker-0.9-cp38-cp38-win_amd64.whl
... ...
invesalius/constants.py
  Diff comments   View file @46e9180
... ... @@ -823,23 +823,35 @@ COIL_ANGLE_ARROW_PROJECTION_THRESHOLD = 5
823 823 CAM_MODE = True
824 824  
825 825 # Tractography visualization
826   -N_TRACTS = 100
827   -PEEL_DEPTH = 5
828   -MAX_PEEL_DEPTH = 30
829   -SEED_OFFSET = 15
  826 +N_TRACTS = 200
  827 +PEEL_DEPTH = 10
  828 +MAX_PEEL_DEPTH = 40
  829 +SEED_OFFSET = 30
830 830 SEED_RADIUS = 1.5
831 831  
832 832 # Increased the default sleep parameter from 0.1 to 0.15 to decrease CPU load during navigation.
833   -SLEEP_NAVIGATION = 0.15
  833 +SLEEP_NAVIGATION = 0.2
834 834 SLEEP_COORDINATES = 0.05
835 835 SLEEP_ROBOT = 0.01
836 836  
837   -BRAIN_OPACITY = 0.5
  837 +BRAIN_OPACITY = 0.6
838 838 N_CPU = psutil.cpu_count()
839   -TREKKER_CONFIG = {'seed_max': 1, 'step_size': 0.1, 'min_fod': 0.1, 'probe_quality': 3,
840   - 'max_interval': 1, 'min_radius_curv': 0.8, 'probe_length': 0.4,
841   - 'write_interval': 50, 'numb_threads': '', 'max_lenth': 200,
842   - 'min_lenth': 20, 'max_sampling_step': 100}
  839 +# the max_sampling_step can be set to something different as well. Above 100 is probably not necessary
  840 +TREKKER_CONFIG = {'seed_max': 1,
  841 + 'step_size': 0.03125,
  842 + 'min_fod': 0.05,
  843 + 'probe_quality': 3,
  844 + 'max_interval': 1,
  845 + 'min_radius_curvature': 0.625,
  846 + 'probe_length': 0.15625,
  847 + 'write_interval': 50,
  848 + 'numb_threads': '',
  849 + 'max_length': 250,
  850 + 'min_length': 10,
  851 + 'max_sampling_step': 100,
  852 + 'data_support_exponent': 0.5,
  853 + 'use_best_init': True,
  854 + 'init_max_est_trials': 100}
843 855  
844 856 MARKER_FILE_MAGICK_STRING = "##INVESALIUS3_MARKER_FILE_"
845 857 CURRENT_MARKER_FILE_VERSION = 0
... ...
invesalius/control.py
  Diff comments   View file @46e9180
... ... @@ -69,7 +69,7 @@ class Controller():
69 69 #DICOM = 1
70 70 #TIFF uCT = 2
71 71 self.img_type = 0
72   - self.affine = None
  72 + self.affine = np.identity(4)
73 73  
74 74 #Init session
75 75 session = ses.Session()
... ... @@ -340,7 +340,7 @@ class Controller():
340 340 if proj.affine:
341 341 self.Slice.affine = np.asarray(proj.affine).reshape(4, 4)
342 342 else:
343   - self.Slice.affine = None
  343 + self.Slice.affine = np.identity(4)
344 344  
345 345 Publisher.sendMessage('Update threshold limits list',
346 346 threshold_range=proj.threshold_range)
... ... @@ -623,6 +623,8 @@ class Controller():
623 623 Publisher.sendMessage(('Set scroll position', 'SAGITAL'),index=proj.matrix_shape[1]/2)
624 624 Publisher.sendMessage(('Set scroll position', 'CORONAL'),index=proj.matrix_shape[2]/2)
625 625  
  626 + # TODO: Check that this is needed with the new way of using affine
  627 + # now the affine should be at least the identity(4) and never None
626 628 if self.Slice.affine is not None:
627 629 Publisher.sendMessage('Enable Go-to-Coord', status=True)
628 630 else:
... ... @@ -731,6 +733,8 @@ class Controller():
731 733 proj.level = self.Slice.window_level
732 734 proj.threshold_range = int(matrix.min()), int(matrix.max())
733 735 proj.spacing = self.Slice.spacing
  736 + # TODO: Check that this is needed with the new way of using affine
  737 + # now the affine should be at least the identity(4) and never None
734 738 if self.Slice.affine is not None:
735 739 proj.affine = self.Slice.affine.tolist()
736 740  
... ...
invesalius/data/coordinates.py
  Diff comments   View file @46e9180
... ... @@ -371,7 +371,7 @@ def DebugCoordRandom(trk_init, trck_id, ref_mode):
371 371 #
372 372 # else:
373 373  
374   - dx = [-70, 70]
  374 + dx = [-30, 30]
375 375 dt = [-180, 180]
376 376  
377 377 coord1 = np.array([uniform(*dx), uniform(*dx), uniform(*dx),
... ...
invesalius/data/coregistration.py
  Diff comments   View file @46e9180
... ... @@ -427,7 +427,7 @@ class CoordinateCorregistrateNoObject(threading.Thread):
427 427 # # seed_aux = pos_world.reshape([1, 4])[0, :3]
428 428 # # seed = seed_aux[np.newaxis, :]
429 429 # #
430   -# # self.tracts = dtr.compute_tracts(tracker, seed, affine_vtk, True)
  430 +# # self.tracts = dtr.compute_and_visualize_tracts(tracker, seed, affine_vtk, True)
431 431 #
432 432 # # wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord)
433 433 # wx.CallAfter(Publisher.sendMessage, 'Update cross position', arg=m_img, position=coord)
... ... @@ -654,7 +654,7 @@ class CoordinateCorregistrateNoObject(threading.Thread):
654 654 # # seed_aux = pos_world.reshape([1, 4])[0, :3]
655 655 # # seed = seed_aux[np.newaxis, :]
656 656 #
657   -# # self.tracts = dtr.compute_tracts(tracker, seed, affine_vtk, True)
  657 +# # self.tracts = dtr.compute_and_visualize_tracts(tracker, seed, affine_vtk, True)
658 658 #
659 659 # # wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord)
660 660 # wx.CallAfter(Publisher.sendMessage, 'Update cross position', arg=m_img, position=coord)
... ...
invesalius/data/styles.py
  Diff comments   View file @46e9180
... ... @@ -544,7 +544,7 @@ class TractsInteractorStyle(CrossInteractorStyle):
544 544  
545 545 def OnTractsMouseClick(self, obj, evt):
546 546 # print("Single mouse click")
547   - # self.tracts = dtr.compute_tracts(self.tracker, self.seed, self.left_pressed)
  547 + # self.tracts = dtr.compute_and_visualize_tracts(self.tracker, self.seed, self.left_pressed)
548 548 self.ChangeTracts(True)
549 549  
550 550 def OnTractsReleaseLeftButton(self, obj, evt):
... ... @@ -554,7 +554,7 @@ class TractsInteractorStyle(CrossInteractorStyle):
554 554  
555 555 def ChangeTracts(self, pressed):
556 556 # print("Trying to compute tracts")
557   - self.tracts = dtr.compute_tracts(self.tracker, self.seed, self.affine_vtk, pressed)
  557 + self.tracts = dtr.compute_and_visualize_tracts(self.tracker, self.seed, self.affine_vtk, pressed)
558 558 # mouse_x, mouse_y = iren.GetEventPosition()
559 559 # wx, wy, wz = self.viewer.get_coordinate_cursor(mouse_x, mouse_y, self.picker)
560 560 # px, py = self.viewer.get_slice_pixel_coord_by_world_pos(wx, wy, wz)
... ...
invesalius/data/surface.py
  Diff comments   View file @46e9180
... ... @@ -365,9 +365,16 @@ class SurfaceManager():
365 365 if self.convert2inv:
366 366 # convert between invesalius and world space with shift in the Y coordinate
367 367 affine = sl.Slice().affine
  368 + # TODO: Check that this is needed with the new way of using affine
  369 + # now the affine should be at least the identity(4) and never None
368 370 if affine is not None:
369   - affine[1, -1] -= sl.Slice().spacing[1] * (sl.Slice().matrix.shape[1] - 1)
  371 + matrix_shape = sl.Slice().matrix.shape
  372 + spacing = sl.Slice().spacing
  373 + img_shift = spacing[1] * (matrix_shape[1] - 1)
  374 + affine = sl.Slice().affine.copy()
  375 + affine[1, -1] -= img_shift
370 376 affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(affine)
  377 +
371 378 actor.SetUserMatrix(affine_vtk)
372 379  
373 380 if overwrite:
... ...
invesalius/data/tractography.py
  Diff comments   View file @46e9180
... ... @@ -29,14 +29,11 @@ import time
29 29 import numpy as np
30 30 import queue
31 31 from invesalius.pubsub import pub as Publisher
32   -from scipy.stats import norm
33 32 import vtk
34 33  
35 34 import invesalius.constants as const
36 35 import invesalius.data.imagedata_utils as img_utils
37 36  
38   -import invesalius.project as prj
39   -
40 37 # Nice print for arrays
41 38 # np.set_printoptions(precision=2)
42 39 # np.set_printoptions(suppress=True)
... ... @@ -104,75 +101,57 @@ def compute_tubes(trk, direction):
104 101 return trk_tube
105 102  
106 103  
107   -def combine_tracts_root(out_list, root, n_block):
  104 +def create_branch(out_list, n_block):
108 105 """Adds a set of tracts to given position in a given vtkMultiBlockDataSet
109 106  
110 107 :param out_list: List of vtkTubeFilters representing the tracts
111 108 :type out_list: list
112   - :param root: A collection of tracts as a vtkMultiBlockDataSet
113   - :type root: vtkMultiBlockDataSet
114 109 :param n_block: The location in the given vtkMultiBlockDataSet to insert the new tracts
115 110 :type n_block: int
116   - :return: The updated collection of tracts as a vtkMultiBlockDataSet
117   - :rtype: vtkMultiBlockDataSet
118   - """
119   -
120   - # create tracts only when at least one was computed
121   - # print("Len outlist in root: ", len(out_list))
122   - if not out_list.count(None) == len(out_list):
123   - for n, tube in enumerate(out_list):
124   - root.SetBlock(n_block + n, tube.GetOutput())
125   -
126   - return root
127   -
128   -
129   -def combine_tracts_branch(out_list):
130   - """Combines a set of tracts in vtkMultiBlockDataSet
131   -
132   - :param out_list: List of vtkTubeFilters representing the tracts
133   - :type out_list: list
134   - :return: A collection of tracts as a vtkMultiBlockDataSet
  111 + :return: The collection of tracts (streamlines) as a vtkMultiBlockDataSet
135 112 :rtype: vtkMultiBlockDataSet
136 113 """
137 114  
  115 + # create a branch and add the streamlines
138 116 branch = vtk.vtkMultiBlockDataSet()
  117 +
139 118 # create tracts only when at least one was computed
140 119 # print("Len outlist in root: ", len(out_list))
  120 + # TODO: check if this if statement is required, because we should
  121 + # call this function only when tracts exist
141 122 if not out_list.count(None) == len(out_list):
142 123 for n, tube in enumerate(out_list):
143   - branch.SetBlock(n, tube.GetOutput())
  124 + branch.SetBlock(n_block + n, tube.GetOutput())
144 125  
145 126 return branch
146 127  
147 128  
148   -def tracts_computation(trk_list, root, n_tracts):
  129 +def compute_tracts(trk_list, n_tract=0, alpha=255):
149 130 """Convert the list of all computed tracts given by Trekker run and returns a vtkMultiBlockDataSet
150 131  
151 132 :param trk_list: List of lists containing the computed tracts and corresponding coordinates
152 133 :type trk_list: list
153   - :param root: A collection of tracts as a vtkMultiBlockDataSet
154   - :type root: vtkMultiBlockDataSet
155   - :param n_tracts:
156   - :type n_tracts: int
  134 + :param n_tract: The integer ID of the block in the vtkMultiBlockDataSet
  135 + :type n_tract: int
  136 + :param alpha: The transparency of the streamlines from 0 to 255 (transparent to opaque)
  137 + :type alpha: int
157 138 :return: The updated collection of tracts as a vtkMultiBlockDataSet
158 139 :rtype: vtkMultiBlockDataSet
159 140 """
160 141  
161 142 # Transform tracts to array
162 143 trk_arr = [np.asarray(trk_n).T if trk_n else None for trk_n in trk_list]
163   -
164 144 # Compute the directions
165   - trk_dir = [compute_directions(trk_n) for trk_n in trk_arr]
166   -
  145 + trk_dir = [compute_directions(trk_n, alpha) for trk_n in trk_arr]
167 146 # Compute the vtk tubes
168 147 out_list = [compute_tubes(trk_arr_n, trk_dir_n) for trk_arr_n, trk_dir_n in zip(trk_arr, trk_dir)]
  148 + # create a branch and add the tracts
  149 + branch = create_branch(out_list, n_tract)
169 150  
170   - root = combine_tracts_root(out_list, root, n_tracts)
171   -
172   - return root
  151 + return branch
173 152  
174 153  
175   -def compute_tracts(trekker, position, affine, affine_vtk, n_tracts):
  154 +def compute_and_visualize_tracts(trekker, position, affine, affine_vtk, n_tracts_max):
176 155 """ Compute tractograms using the Trekker library.
177 156  
178 157 :param trekker: Trekker library instance
... ... @@ -183,52 +162,51 @@ def compute_tracts(trekker, position, affine, affine_vtk, n_tracts):
183 162 :type affine: numpy.ndarray
184 163 :param affine_vtk: vtkMatrix4x4 isntance with affine transformation matrix
185 164 :type affine_vtk: vtkMatrix4x4
186   - :param n_tracts: number of tracts to compute
187   - :type n_tracts: int
  165 + :param n_tracts_max: maximum number of tracts to compute
  166 + :type n_tracts_max: int
188 167 """
189 168  
190   - # during neuronavigation, root needs to be initialized outside the while loop so the new tracts
191   - # can be appended to the root block set
192   - root = vtk.vtkMultiBlockDataSet()
  169 + # root = vtk.vtkMultiBlockDataSet()
193 170 # Juuso's
194 171 # seed = np.array([[-8.49, -8.39, 2.5]])
195 172 # Baran M1
196 173 # seed = np.array([[27.53, -77.37, 46.42]])
197 174 seed_trk = img_utils.convert_world_to_voxel(position, affine)
198   - # print("seed example: {}".format(seed_trk))
199   - trekker.seed_coordinates(np.repeat(seed_trk, n_tracts, axis=0))
200   - # print("trk list len: ", len(trekker.run()))
201   - trk_list = trekker.run()
202   - if trk_list:
203   - root = tracts_computation(trk_list, root, 0)
204   - Publisher.sendMessage('Remove tracts')
205   - Publisher.sendMessage('Update tracts', root=root, affine_vtk=affine_vtk, coord_offset=position)
206   - else:
207   - Publisher.sendMessage('Remove tracts')
208   -
209   -
210   -def tracts_computation_branch(trk_list, alpha=255):
211   - """Convert the list of all computed tracts given by Trekker run and returns a vtkMultiBlockDataSet
212   -
213   - :param trk_list: List of lists containing the computed tracts and corresponding coordinates
214   - :type trk_list: list
215   - :param alpha: opacity value in the interval [0, 255]. The 0 is no opacity (total transparency).
216   - :type alpha: int
217   - :return: The collection of tracts as a vtkMultiBlockDataSet
218   - :rtype: vtkMultiBlockDataSet
219   - """
220   - # Transform tracts to array
221   - trk_arr = [np.asarray(trk_n).T if trk_n else None for trk_n in trk_list]
222   - # Compute the directions
223   - trk_dir = [compute_directions(trk_n, alpha) for trk_n in trk_arr]
224   - # Compute the vtk tubes
225   - tube_list = [compute_tubes(trk_arr_n, trk_dir_n) for trk_arr_n, trk_dir_n in zip(trk_arr, trk_dir)]
226   - branch = combine_tracts_branch(tube_list)
227   -
228   - return branch
  175 + bundle = vtk.vtkMultiBlockDataSet()
  176 + n_branches, n_tracts, count_loop = 0, 0, 0
  177 + n_threads = 2 * const.N_CPU - 1
  178 +
  179 + while n_tracts < n_tracts_max:
  180 + n_param = 1 + (count_loop % 10)
  181 + # rescale the alpha value that defines the opacity of the branch
  182 + # the n interval is [1, 10] and the new interval is [51, 255]
  183 + # the new interval is defined to have no 0 opacity (minimum is 51, i.e., 20%)
  184 + alpha = (n_param - 1) * (255 - 51) / (10 - 1) + 51
  185 + trekker.minFODamp(n_param * 0.01)
  186 +
  187 + # print("seed example: {}".format(seed_trk))
  188 + trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0))
  189 + # print("trk list len: ", len(trekker.run()))
  190 + trk_list = trekker.run()
  191 + n_tracts += len(trk_list)
  192 + if len(trk_list):
  193 + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha)
  194 + bundle.SetBlock(n_branches, branch)
  195 + n_branches += 1
  196 +
  197 + count_loop += 1
  198 +
  199 + if (count_loop == 20) and (n_tracts == 0):
  200 + break
  201 +
  202 + Publisher.sendMessage('Remove tracts')
  203 + if n_tracts:
  204 + Publisher.sendMessage('Update tracts', root=bundle, affine_vtk=affine_vtk,
  205 + coord_offset=position, coord_offset_w=seed_trk[0].tolist())
229 206  
230 207  
231 208 class ComputeTractsThread(threading.Thread):
  209 + # TODO: Remove this class and create a case where no ACT is provided in the class ComputeTractsACTThread
232 210  
233 211 def __init__(self, inp, queues, event, sle):
234 212 """Class (threading) to compute real time tractography data for visualization.
... ... @@ -265,6 +243,7 @@ class ComputeTractsThread(threading.Thread):
265 243  
266 244 trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.inp
267 245 # n_threads = n_tracts_total
  246 + n_threads = int(n_threads/4)
268 247 p_old = np.array([[0., 0., 0.]])
269 248 n_tracts = 0
270 249  
... ... @@ -312,13 +291,13 @@ class ComputeTractsThread(threading.Thread):
312 291 # run the trekker, this is the slowest line of code, be careful to just use once!
313 292 trk_list = trekker.run()
314 293  
315   - if trk_list:
  294 + if len(trk_list) > 2:
316 295 # print("dist: {}".format(dist))
317 296 if dist >= seed_radius:
318 297 # when moving the coil further than the seed_radius restart the bundle computation
319 298 bundle = vtk.vtkMultiBlockDataSet()
320 299 n_branches = 0
321   - branch = tracts_computation_branch(trk_list)
  300 + branch = compute_tracts(trk_list, n_tract=0, alpha=255)
322 301 bundle.SetBlock(n_branches, branch)
323 302 n_branches += 1
324 303 n_tracts = branch.GetNumberOfBlocks()
... ... @@ -326,7 +305,7 @@ class ComputeTractsThread(threading.Thread):
326 305 # TODO: maybe keep computing even if reaches the maximum
327 306 elif dist < seed_radius and n_tracts < n_tracts_total:
328 307 # compute tracts blocks and add to bungle until reaches the maximum number of tracts
329   - branch = tracts_computation_branch(trk_list)
  308 + branch = compute_tracts(trk_list, n_tract=0, alpha=255)
330 309 if bundle:
331 310 bundle.SetBlock(n_branches, branch)
332 311 n_tracts += branch.GetNumberOfBlocks()
... ... @@ -361,314 +340,218 @@ class ComputeTractsThread(threading.Thread):
361 340  
362 341 class ComputeTractsACTThread(threading.Thread):
363 342  
364   - def __init__(self, inp, queues, event, sle):
  343 + def __init__(self, input_list, queues, event, sleep_thread):
365 344 """Class (threading) to compute real time tractography data for visualization.
366 345  
367 346 Tracts are computed using the Trekker library by Baran Aydogan (https://dmritrekker.github.io/)
368 347 For VTK visualization, each tract (fiber) is a constructed as a tube and many tubes combined in one
369 348 vtkMultiBlockDataSet named as a branch. Several branches are combined in another vtkMultiBlockDataSet named as
370 349 bundle, to obtain fast computation and visualization. The bundle dataset is mapped to a single vtkActor.
371   - Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the invesalius 3D scene.
  350 + Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the
  351 + invesalius 3D scene.
372 352  
373 353 Sleep function in run method is used to avoid blocking GUI and more fluent, real-time navigation
374 354  
375   - :param inp: List of inputs: trekker instance, affine numpy array, seed_offset, seed_radius, n_threads
376   - :type inp: list
  355 + :param input_list: List of inputs: trekker instance, affine numpy array, seed offset, total number of tracts,
  356 + seed radius, number of threads in computer, ACT data array, affine vtk matrix,
  357 + image shift for vtk to mri transformation
  358 + :type input_list: list
377 359 :param queues: Queue list with coord_tracts_queue (Queue instance that manage co-registered coordinates) and
378 360 tracts_queue (Queue instance that manage the tracts to be visualized)
379 361 :type queues: list[queue.Queue, queue.Queue]
380 362 :param event: Threading event to coordinate when tasks as done and allow UI release
381 363 :type event: threading.Event
382   - :param sle: Sleep pause in seconds
383   - :type sle: float
  364 + :param sleep_thread: Sleep pause in seconds
  365 + :type sleep_thread: float
384 366 """
385 367  
386 368 threading.Thread.__init__(self, name='ComputeTractsThreadACT')
387   - self.inp = inp
388   - # self.coord_queue = coord_queue
  369 + self.input_list = input_list
389 370 self.coord_tracts_queue = queues[0]
390 371 self.tracts_queue = queues[1]
391   - # on first pilots (january 12, 2021) used (-4, 4)
392   - self.coord_list_w = img_utils.create_grid((-2, 2), (0, 20), inp[2]-5, 1)
393   - # self.coord_list_sph = img_utils.create_spherical_grid(10, 1)
394   - # self.coord_list_sph = img_utils.create_spherical_grid(10, 1)
395   - # x_norm = np.linspace(norm.ppf(0.01), norm.ppf(0.99), 2*self.coord_list_sph.shape[0])
396   - # self.pdf = np.flipud(norm.pdf(x_norm[:self.coord_list_sph.shape[0]], loc=0, scale=2.))
397   - # self.sph_idx = np.linspace(0, self.coord_list_sph.shape[0], num=self.coord_list_sph.shape[0], dtype=int)
398   - # self.visualization_queue = visualization_queue
399 372 self.event = event
400   - self.sle = sle
401   -
402   - # prj_data = prj.Project()
403   - # matrix_shape = tuple(prj_data.matrix_shape)
404   - # self.img_shift = matrix_shape[1]
  373 + self.sleep_thread = sleep_thread
405 374  
406 375 def run(self):
407 376  
408   - # trekker, affine, offset, n_tracts_total, seed_radius, n_threads = self.inp
409   - trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.inp
  377 + trekker, affine, offset, n_tracts_total, seed_radius, n_threads, act_data, affine_vtk, img_shift = self.input_list
410 378  
411   - # n_threads = n_tracts_total
412 379 p_old = np.array([[0., 0., 0.]])
413   - p_old_pre = np.array([[0., 0., 0.]])
414   - coord_offset = None
415   - n_tracts = 0
416   - n_branches = 0
  380 + n_branches, n_tracts, count_loop = 0, 0, 0
417 381 bundle = None
418   - sph_sampling = True
419 382 dist_radius = 1.5
420 383  
421   - xyz = img_utils.random_sample_sphere(radius=seed_radius, size=100)
422   - coord_list_sph = np.hstack([xyz, np.ones([xyz.shape[0], 1])]).T
  384 + # TODO: Try a denser and bigger grid, because it's just a matrix multiplication
  385 + # maybe 15 mm below the coil offset by default and 5 cm deep
  386 + # create the rectangular grid to find the gray-white matter boundary
  387 + coord_list_w = img_utils.create_grid((-2, 2), (0, 20), offset - 5, 1)
  388 +
  389 + # create the spherical grid to sample around the seed location
  390 + samples_in_sphere = img_utils.random_sample_sphere(radius=seed_radius, size=100)
  391 + coord_list_sphere = np.hstack([samples_in_sphere, np.ones([samples_in_sphere.shape[0], 1])]).T
423 392 m_seed = np.identity(4)
  393 +
424 394 # Compute the tracts
425   - # print('ComputeTractsThread: event {}'.format(self.event.is_set()))
426 395 while not self.event.is_set():
427 396 try:
428   - # print("Computing tracts")
429 397 # get from the queue the coordinates, coregistration transformation matrix, and flipped matrix
430 398 m_img_flip = self.coord_tracts_queue.get_nowait()
431   - # coord, m_img, m_img_flip = self.coord_queue.get_nowait()
432   - # print('ComputeTractsThread: get {}'.format(count))
433   -
434   - # TODO: Remove this is not needed
435   - # 20200402: in this new refactored version the m_img comes different than the position
436   - # the new version m_img is already flixped in y, which means that Y is negative
437   - # if only the Y is negative maybe no need for the flip_x funtcion at all in the navigation
438   - # but check all coord_queue before why now the m_img comes different than position
439   - # 20200403: indeed flip_x is just a -1 multiplication to the Y coordinate, remove function flip_x
440   - # m_img_flip = m_img.copy()
441   - # m_img_flip[1, -1] = -m_img_flip[1, -1]
442 399  
443 400 # DEBUG: Uncomment the m_img_flip below so that distance is fixed and tracts keep computing
444 401 # m_img_flip[:3, -1] = (5., 10., 12.)
445   - dist = abs(np.linalg.norm(p_old_pre - np.asarray(m_img_flip[:3, -1])))
446   - p_old_pre = m_img_flip[:3, -1].copy()
447   -
448   - # Uncertanity visualization --
449   - # each tract branch is computed with one set of parameters ajusted from 1 to 10
450   - n_param = 1 + (n_branches % 10)
  402 + dist = abs(np.linalg.norm(p_old - np.asarray(m_img_flip[:3, -1])))
  403 + p_old = m_img_flip[:3, -1].copy()
  404 +
  405 + # Uncertainty visualization --
  406 + # each tract branch is computed with one minFODamp adjusted from 0.01 to 0.1
  407 + # the higher the minFODamp the more the streamlines are faithful to the data, so they become more strict
  408 + # but also may loose some connections.
  409 + # the lower the more relaxed streamline also increases the chance of false positives
  410 + n_param = 1 + (count_loop % 10)
451 411 # rescale the alpha value that defines the opacity of the branch
452 412 # the n interval is [1, 10] and the new interval is [51, 255]
453 413 # the new interval is defined to have no 0 opacity (minimum is 51, i.e., 20%)
454 414 alpha = (n_param - 1) * (255 - 51) / (10 - 1) + 51
455 415 trekker.minFODamp(n_param * 0.01)
456   - trekker.dataSupportExponent(n_param * 0.1)
457 416 # ---
458 417  
  418 + try:
  419 + # The original seed location is replaced by the gray-white matter interface that is closest to
  420 + # the coil center
  421 + coord_list_w_tr = m_img_flip @ coord_list_w
  422 + coord_offset = grid_offset(act_data, coord_list_w_tr, img_shift)
  423 + except IndexError:
  424 + # This error might be caused by the coordinate exceeding the image array dimensions.
  425 + # This happens during navigation because the coil location can have coordinates outside the image
  426 + # boundaries
  427 + # Translate the coordinate along the normal vector of the object/coil
  428 + coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2]
  429 + # ---
  430 +
  431 + # Spherical sampling of seed coordinates ---
  432 + # compute the samples of a sphere centered on seed coordinate offset by the grid
  433 + # given in the invesalius-vtk space
  434 + samples = np.random.choice(coord_list_sphere.shape[1], size=100)
  435 + m_seed[:-1, -1] = coord_offset.copy()
  436 + # translate the spherical grid samples to the coil location in invesalius-vtk space
  437 + seed_trk_r_inv = m_seed @ coord_list_sphere[:, samples]
  438 +
  439 + coord_offset_w = np.hstack((coord_offset, 1.0)).reshape([4, 1])
  440 +
  441 + try:
  442 + # Anatomic constrained seed computation ---
  443 + # find only the samples inside the white matter as with the ACT enabled in the Trekker,
  444 + # it does not compute any tracts outside the white matter
  445 + # convert from inveslaius-vtk to mri space
  446 + seed_trk_r_mri = seed_trk_r_inv[:3, :].T.astype(int) + np.array([[0, img_shift, 0]], dtype=np.int32)
  447 + labs = act_data[seed_trk_r_mri[..., 0], seed_trk_r_mri[..., 1], seed_trk_r_mri[..., 2]]
  448 + # find all samples in the white matter
  449 + labs_id = np.where(labs == 1)
  450 + # Trekker has internal multiprocessing approach done in C. Here the number of available threads - 1
  451 + # is given, but in case a large number of tracts is requested, it will compute all in parallel
  452 + # automatically for a more fluent navigation, better to compute the maximum number the computer can
  453 + # handle otherwise gets too slow for the multithreading in Python
  454 + seed_trk_r_inv_sampled = seed_trk_r_inv[:, labs_id[0][:n_threads]]
  455 +
  456 + except IndexError:
  457 + # same as on the grid offset above, if the coil is too far from the mri volume the array indices
  458 + # are beyond the mri boundaries
  459 + # in this case use the grid center instead of the spherical samples
  460 + seed_trk_r_inv_sampled = coord_offset_w.copy()
  461 +
  462 + # convert to the world coordinate system for trekker
  463 + seed_trk_r_world_sampled = np.linalg.inv(affine) @ seed_trk_r_inv_sampled
  464 + seed_trk_r_world_sampled = seed_trk_r_world_sampled.T[:, :3]
  465 +
  466 + # convert to the world coordinate system for saving in the marker list
  467 + coord_offset_w = np.linalg.inv(affine) @ coord_offset_w
  468 + coord_offset_w = np.squeeze(coord_offset_w.T[:, :3])
  469 +
  470 + # DEBUG: uncomment the seed_trk below
  471 + # seed_trk.shape == [1, 3]
  472 + # Juuso's
  473 + # seed_trk = np.array([[-8.49, -8.39, 2.5]])
  474 + # Baran M1
  475 + # seed_trk = np.array([[27.53, -77.37, 46.42]])
  476 + # print("Given: {}".format(seed_trk.shape))
  477 + # print("Seed: {}".format(seed))
  478 + # joonas seed that has good tracts
  479 + # seed_trk = np.array([[29.12, -13.33, 31.65]])
  480 + # seed_trk_img = np.array([[117, 127, 161]])
  481 +
459 482 # When moving the coil further than the seed_radius restart the bundle computation
460 483 # Currently, it stops to compute tracts when the maximum number of tracts is reached maybe keep
461 484 # computing even if reaches the maximum
462 485 if dist >= dist_radius:
463   - # Anatomic constrained seed computation ---
464   - # The original seed location is replaced by the gray-white matter interface that is closest to
465   - # the coil center
466   - try:
467   - #TODO: Create a dialog error to say when the ACT data is not loaded and prevent
468   - # the interface from freezing. Give the user a chance to load it (maybe in task_navigator)
469   - coord_list_w_tr = m_img_flip @ self.coord_list_w
470   - coord_offset = grid_offset(act_data, coord_list_w_tr, img_shift)
471   - except IndexError:
472   - # This error might be caused by the coordinate exceeding the image array dimensions.
473   - # Needs further verification.
474   - coord_offset = None
475   - # ---
476   -
477   - # Translate the coordinate along the normal vector of the object/coil ---
478   - if coord_offset is None:
479   - # apply the coil transformation matrix
480   - coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2]
481   - # ---
482   -
483   - # convert the world coordinates to the voxel space for using as a seed in Trekker
484   - # seed_trk.shape == [1, 3]
485   - seed_trk = img_utils.convert_world_to_voxel(coord_offset, affine)
486   - # print("Desired: {}".format(seed_trk.shape))
487   -
488   - # DEBUG: uncomment the seed_trk below
489   - # Juuso's
490   - # seed_trk = np.array([[-8.49, -8.39, 2.5]])
491   - # Baran M1
492   - # seed_trk = np.array([[27.53, -77.37, 46.42]])
493   - # print("Given: {}".format(seed_trk.shape))
494   - # print("Seed: {}".format(seed))
495   -
496   - # Spherical sampling of seed coordinates ---
497   - if sph_sampling:
498   - # CHECK: We use ACT only for the origin seed, but not for all the other coordinates.
499   - # Check how this can be solved. Applying ACT to all coordinates is maybe too much.
500   - # Maybe it doesn't matter because the ACT is just to help finding the closest location to
501   - # the TMS coil center. Also, note that the spherical sampling is applied only when the coil
502   - # location changes, all further iterations used the fixed n_threads samples to compute the
503   - # remaining tracts.
504   -
505   - # samples = np.random.choice(self.sph_idx, size=n_threads, p=self.pdf)
506   - # m_seed[:-1, -1] = seed_trk
507   - # sph_seed = m_seed @ self.coord_list_sph
508   - # seed_trk_r = sph_seed[samples, :]
509   - samples = np.random.choice(coord_list_sph.shape[1], size=n_threads)
510   - m_seed[:-1, -1] = seed_trk
511   - seed_trk_r = m_seed @ coord_list_sph[:, samples]
512   - seed_trk_r = seed_trk_r[:-1, :].T
513   - else:
514   - # set the seeds for trekker, one seed is repeated n_threads times
515   - # shape (24, 3)
516   - seed_trk_r = np.repeat(seed_trk, n_threads, axis=0)
517   -
518   - # ---
519   -
520   - # Trekker has internal multiprocessing approach done in C. Here the number of available threads is
521   - # given, but in case a large number of tracts is requested, it will compute all in parallel
522   - # automatically for a more fluent navigation, better to compute the maximum number the computer
523   - # handles
524   - trekker.seed_coordinates(seed_trk_r)
525   - # trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0))
  486 + bundle = None
  487 + n_tracts, n_branches = 0, 0
526 488  
  489 + # we noticed that usually the navigation lags or crashes when moving the coil location
  490 + # to reduce the overhead for when the coil is moving, we compute only half the number of tracts
  491 + # that should be computed when the coil is fixed in the same location
  492 + # required input is Nx3 array
  493 + trekker.seed_coordinates(seed_trk_r_world_sampled[::2, :])
527 494 # run the trekker, this is the slowest line of code, be careful to just use once!
528 495 trk_list = trekker.run()
529 496  
530 497 # check if any tract was found, otherwise doesn't count
531   - if trk_list:
  498 + if len(trk_list):
  499 + # a bundle consists for multiple branches and each branch consists of multiple streamlines
  500 + # every iteration in the main loop adds a branch to the bundle
  501 + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha)
  502 + n_tracts = branch.GetNumberOfBlocks()
  503 +
  504 + # create and add branch to the bundle
532 505 bundle = vtk.vtkMultiBlockDataSet()
533   - branch = tracts_computation_branch(trk_list, alpha)
534 506 bundle.SetBlock(n_branches, branch)
535 507 n_branches = 1
536   - n_tracts = branch.GetNumberOfBlocks()
537   - else:
538   - bundle = None
539   - n_branches = 0
540   - n_tracts = 0
541 508  
542 509 elif dist < dist_radius and n_tracts < n_tracts_total:
  510 + # when the coil is fixed in place and the number of tracts is smaller than the total
  511 + if not bundle:
  512 + # same as above, when creating the bundle (vtkMultiBlockDataSet) we only compute half the number
  513 + # of tracts to reduce the overhead
  514 + bundle = vtk.vtkMultiBlockDataSet()
  515 + # required input is Nx3 array
  516 + trekker.seed_coordinates(seed_trk_r_world_sampled[::2, :])
  517 + n_tracts, n_branches = 0, 0
  518 + else:
  519 + # if the bundle exists compute all tracts requested
  520 + # required input is Nx3 array
  521 + trekker.seed_coordinates(seed_trk_r_world_sampled)
  522 +
543 523 trk_list = trekker.run()
544   - if trk_list:
  524 +
  525 + if len(trk_list):
545 526 # compute tract blocks and add to bundle until reaches the maximum number of tracts
546 527 # the alpha changes depending on the parameter set
547   - branch = tracts_computation_branch(trk_list, alpha)
548   - if bundle:
549   - bundle.SetBlock(n_branches, branch)
550   - n_tracts += branch.GetNumberOfBlocks()
551   - n_branches += 1
552   - else:
553   - bundle = vtk.vtkMultiBlockDataSet()
554   - bundle.SetBlock(n_branches, branch)
555   - n_branches = 1
556   - n_tracts = branch.GetNumberOfBlocks()
557   - # else:
558   - # bundle = None
  528 + branch = compute_tracts(trk_list, n_tract=0, alpha=alpha)
  529 + n_tracts += branch.GetNumberOfBlocks()
  530 + # add branch to the bundle
  531 + bundle.SetBlock(n_branches, branch)
  532 + n_branches += 1
559 533  
560   - # else:
561   - # bundle = None
  534 + # keep adding to the number of loops even if the tracts were not find
  535 + # this will keep the minFODamp changing and new seed coordinates being tried which would allow
  536 + # higher probability of finding tracts
  537 + count_loop += 1
562 538  
563   - # rethink if this should be inside the if condition, it may lock the thread if no tracts are found
564   - # use no wait to ensure maximum speed and avoid visualizing old tracts in the queue, this might
  539 + # use "nowait" to ensure maximum speed and avoid visualizing old tracts in the queue, this might
565 540 # be more evident in slow computer or for heavier tract computations, it is better slow update
566 541 # than visualizing old data
567   - # self.visualization_queue.put_nowait([coord, m_img, bundle])
568   - self.tracts_queue.put_nowait((bundle, affine_vtk, coord_offset))
569   - # print('ComputeTractsThread: put {}'.format(count))
570   -
  542 + self.tracts_queue.put_nowait((bundle, affine_vtk, coord_offset, coord_offset_w))
571 543 self.coord_tracts_queue.task_done()
572   - # self.coord_queue.task_done()
573   - # print('ComputeTractsThread: done {}'.format(count))
574 544  
575 545 # sleep required to prevent user interface from being unresponsive
576   - time.sleep(self.sle)
  546 + time.sleep(self.sleep_thread)
577 547 # if no coordinates pass
578 548 except queue.Empty:
579   - # print("Empty queue in tractography")
580 549 pass
581 550 # if queue is full mark as done (may not be needed in this new "nowait" method)
582 551 except queue.Full:
583   - # self.coord_queue.task_done()
584 552 self.coord_tracts_queue.task_done()
585 553  
586 554  
587   -class ComputeTractsThreadSingleBlock(threading.Thread):
588   -
589   - def __init__(self, inp, affine_vtk, coord_queue, visualization_queue, event, sle):
590   - """Class (threading) to compute real time tractography data for visualization in a single loop.
591   -
592   - Different than ComputeTractsThread because it does not keep adding tracts to the bundle until maximum,
593   - is reached. It actually compute all requested tracts at once. (Might be deleted in the future)!
594   - Tracts are computed using the Trekker library by Baran Aydogan (https://dmritrekker.github.io/)
595   - For VTK visualization, each tract (fiber) is a constructed as a tube and many tubes combined in one
596   - vtkMultiBlockDataSet named as a branch. Several branches are combined in another vtkMultiBlockDataSet named as
597   - bundle, to obtain fast computation and visualization. The bundle dataset is mapped to a single vtkActor.
598   - Mapper and Actor are computer in the data/viewer_volume.py module for easier handling in the invesalius 3D scene.
599   -
600   - Sleep function in run method is used to avoid blocking GUI and more fluent, real-time navigation
601   -
602   - :param inp: List of inputs: trekker instance, affine numpy array, seed_offset, seed_radius, n_threads
603   - :type inp: list
604   - :param affine_vtk: Affine matrix in vtkMatrix4x4 instance to update objects position in 3D scene
605   - :type affine_vtk: vtkMatrix4x4
606   - :param coord_queue: Queue instance that manage coordinates read from tracking device and coregistered
607   - :type coord_queue: queue.Queue
608   - :param visualization_queue: Queue instance that manage coordinates to be visualized
609   - :type visualization_queue: queue.Queue
610   - :param event: Threading event to coordinate when tasks as done and allow UI release
611   - :type event: threading.Event
612   - :param sle: Sleep pause in seconds
613   - :type sle: float
614   - """
615   -
616   - threading.Thread.__init__(self, name='ComputeTractsThread')
617   - self.inp = inp
618   - self.affine_vtk = affine_vtk
619   - self.coord_queue = coord_queue
620   - self.visualization_queue = visualization_queue
621   - self.event = event
622   - self.sle = sle
623   -
624   - def run(self):
625   -
626   - trekker, affine, offset, n_tracts_total, seed_radius, n_threads = self.inp
627   - # as a single block, computes all the maximum number of tracts at once, not optimal for navigation
628   - n_threads = n_tracts_total
629   - p_old = np.array([[0., 0., 0.]])
630   - root = vtk.vtkMultiBlockDataSet()
631   -
632   - # Compute the tracts
633   - # print('ComputeTractsThread: event {}'.format(self.event.is_set()))
634   - while not self.event.is_set():
635   - try:
636   - coord, [coord_raw, markers_flag], m_img, m_img_flip = self.coord_queue.get_nowait()
637   -
638   - # translate the coordinate along the normal vector of the object/coil
639   - coord_offset = m_img_flip[:3, -1] - offset * m_img_flip[:3, 2]
640   - # coord_offset = np.array([[27.53, -77.37, 46.42]])
641   - dist = abs(np.linalg.norm(p_old - np.asarray(coord_offset)))
642   - p_old = coord_offset.copy()
643   - seed_trk = img_utils.convert_world_to_voxel(coord_offset, affine)
644   - # Juuso's
645   - # seed_trk = np.array([[-8.49, -8.39, 2.5]])
646   - # Baran M1
647   - # seed_trk = np.array([[27.53, -77.37, 46.42]])
648   - # print("Seed: {}".format(seed))
649   -
650   - # set the seeds for trekker, one seed is repeated n_threads times
651   - # trekker has internal multiprocessing approach done in C. Here the number of available threads is give,
652   - # but in case a large number of tracts is requested, it will compute all in parallel automatically
653   - # for a more fluent navigation, better to compute the maximum number the computer handles
654   - trekker.seed_coordinates(np.repeat(seed_trk, n_threads, axis=0))
655   - # run the trekker, this is the slowest line of code, be careful to just use once!
656   - trk_list = trekker.run()
657   -
658   - if trk_list:
659   - # if the seed is outside the defined radius, restart the bundle computation
660   - if dist >= seed_radius:
661   - root = tracts_computation(trk_list, root, 0)
662   - self.visualization_queue.put_nowait((coord, m_img, root))
663   -
664   - self.coord_queue.task_done()
665   - time.sleep(self.sle)
666   - except queue.Empty:
667   - pass
668   - except queue.Full:
669   - self.coord_queue.task_done()
670   -
671   -
672 555 def set_trekker_parameters(trekker, params):
673 556 """Set all user-defined parameters for tractography computation using the Trekker library
674 557  
... ... @@ -680,21 +563,27 @@ def set_trekker_parameters(trekker, params):
680 563 :rtype: list
681 564 """
682 565 trekker.seed_maxTrials(params['seed_max'])
683   - # trekker.stepSize(params['step_size'])
684   - trekker.minFODamp(params['min_fod'])
685   - # trekker.probeQuality(params['probe_quality'])
686   - # trekker.maxEstInterval(params['max_interval'])
687   - # trekker.minRadiusOfCurvature(params['min_radius_curv'])
688   - # trekker.probeLength(params['probe_length'])
  566 + trekker.stepSize(params['step_size'])
  567 + # minFODamp is not set because it should vary in the loop to create the
  568 + # different transparency tracts
  569 + # trekker.minFODamp(params['min_fod'])
  570 + trekker.probeQuality(params['probe_quality'])
  571 + trekker.maxEstInterval(params['max_interval'])
  572 + trekker.minRadiusOfCurvature(params['min_radius_curvature'])
  573 + trekker.probeLength(params['probe_length'])
689 574 trekker.writeInterval(params['write_interval'])
690   - trekker.maxLength(params['max_lenth'])
691   - trekker.minLength(params['min_lenth'])
  575 + # these two does not need to be set in the new package
  576 + # trekker.maxLength(params['max_length'])
  577 + trekker.minLength(params['min_length'])
692 578 trekker.maxSamplingPerStep(params['max_sampling_step'])
  579 + trekker.dataSupportExponent(params['data_support_exponent'])
  580 + #trekker.useBestAtInit(params['use_best_init'])
  581 + #trekker.initMaxEstTrials(params['init_max_est_trials'])
693 582  
694 583 # check number if number of cores is valid in configuration file,
695 584 # otherwise use the maximum number of threads which is usually 2*N_CPUS
696   - n_threads = 2 * const.N_CPU
697   - if isinstance((params['numb_threads']), int) and params['numb_threads'] <= 2*const.N_CPU:
  585 + n_threads = 2 * const.N_CPU - 1
  586 + if isinstance((params['numb_threads']), int) and params['numb_threads'] <= (2*const.N_CPU-1):
698 587 n_threads = params['numb_threads']
699 588  
700 589 trekker.numberOfThreads(n_threads)
... ... @@ -712,13 +601,21 @@ def grid_offset(data, coord_list_w_tr, img_shift):
712 601  
713 602 # extract the first occurrence of a specific label from the MRI image
714 603 labs = data[coord_list_w_tr_mri[..., 0], coord_list_w_tr_mri[..., 1], coord_list_w_tr_mri[..., 2]]
715   - lab_first = np.argmax(labs == 1)
716   - if labs[lab_first] == 1:
717   - pt_found = coord_list_w_tr_mri[lab_first, :]
  604 + lab_first = np.where(labs == 1)
  605 + if not lab_first:
  606 + pt_found_inv = None
  607 + else:
  608 + pt_found = coord_list_w_tr[:, lab_first[0][0]][:3]
718 609 # convert coordinate back to invesalius 3D space
719 610 pt_found_inv = pt_found - np.array([0., img_shift, 0.])
720   - else:
721   - pt_found_inv = None
  611 +
  612 + # lab_first = np.argmax(labs == 1)
  613 + # if labs[lab_first] == 1:
  614 + # pt_found = coord_list_w_tr_mri[lab_first, :]
  615 + # # convert coordinate back to invesalius 3D space
  616 + # pt_found_inv = pt_found - np.array([0., img_shift, 0.])
  617 + # else:
  618 + # pt_found_inv = None
722 619  
723 620 # # convert to world coordinate space to use as seed for fiber tracking
724 621 # pt_found_tr = np.append(pt_found, 1)[np.newaxis, :].T
... ...
invesalius/data/viewer_volume.py
  Diff comments   View file @46e9180
... ... @@ -307,7 +307,6 @@ class Viewer(wx.Panel):
307 307 Publisher.subscribe(self.OnRemoveTracts, 'Remove tracts')
308 308 Publisher.subscribe(self.UpdateSeedOffset, 'Update seed offset')
309 309 Publisher.subscribe(self.UpdateMarkerOffsetState, 'Update marker offset state')
310   - Publisher.subscribe(self.UpdateMarkerOffsetPosition, 'Update marker offset')
311 310 Publisher.subscribe(self.AddPeeledSurface, 'Update peel')
312 311 Publisher.subscribe(self.GetPeelCenters, 'Get peel centers and normals')
313 312 Publisher.subscribe(self.Initlocator_viewer, 'Get init locator')
... ... @@ -857,8 +856,10 @@ class Viewer(wx.Panel):
857 856 else:
858 857 self.DisableCoilTracker()
859 858 if self.actor_peel:
860   - self.object_orientation_torus_actor.SetVisibility(1)
861   - self.obj_projection_arrow_actor.SetVisibility(1)
  859 + if self.object_orientation_torus_actor:
  860 + self.object_orientation_torus_actor.SetVisibility(1)
  861 + if self.obj_projection_arrow_actor:
  862 + self.obj_projection_arrow_actor.SetVisibility(1)
862 863  
863 864 def OnUpdateObjectTargetGuide(self, m_img, coord):
864 865  
... ... @@ -1607,10 +1608,6 @@ class Viewer(wx.Panel):
1607 1608 self.ren.AddActor(self.mark_actor)
1608 1609 self.Refresh()
1609 1610  
1610   - def UpdateMarkerOffsetPosition(self, coord_offset):
1611   - self.mark_actor.SetPosition(coord_offset)
1612   - self.Refresh()
1613   -
1614 1611 def UpdateObjectOrientation(self, m_img, coord):
1615 1612 # print("Update object orientation")
1616 1613  
... ... @@ -1687,7 +1684,7 @@ class Viewer(wx.Panel):
1687 1684 # self.ball_actor.SetVisibility(1)
1688 1685 self.Refresh()
1689 1686  
1690   - def OnUpdateTracts(self, root=None, affine_vtk=None, coord_offset=None):
  1687 + def OnUpdateTracts(self, root=None, affine_vtk=None, coord_offset=None, coord_offset_w=None):
1691 1688 mapper = vtk.vtkCompositePolyDataMapper2()
1692 1689 mapper.SetInputDataObject(root)
1693 1690  
... ...
invesalius/gui/dialogs.py
  Diff comments   View file @46e9180
... ... @@ -4203,7 +4203,7 @@ class GoToDialogScannerCoord(wx.Dialog):
4203 4203 main_sizer.Fit(self)
4204 4204  
4205 4205 self.orientation = None
4206   - self.affine = None
  4206 + self.affine = np.identity(4)
4207 4207  
4208 4208 self.__bind_events()
4209 4209  
... ...
invesalius/gui/task_navigator.py
  Diff comments   View file @46e9180
... ... @@ -165,7 +165,7 @@ class InnerFoldPanel(wx.Panel):
165 165 # Study this.
166 166  
167 167 fold_panel = fpb.FoldPanelBar(self, -1, wx.DefaultPosition,
168   - (10, 310), 0, fpb.FPB_SINGLE_FOLD)
  168 + (10, 330), 0, fpb.FPB_SINGLE_FOLD)
169 169  
170 170 # Initialize Tracker and PedalConnection objects here so that they are available to several panels.
171 171 #
... ... @@ -1078,17 +1078,15 @@ class ObjectRegistrationPanel(wx.Panel):
1078 1078  
1079 1079 try:
1080 1080 if filename:
1081   - #TODO: Improve method to read the file, using "with" similar to OnLoadParameters
1082   - data = np.loadtxt(filename, delimiter='\t')
1083   - self.obj_fiducials = data[:, :3]
1084   - self.obj_orients = data[:, 3:]
  1081 + with open(filename, 'r') as text_file:
  1082 + data = [s.split('\t') for s in text_file.readlines()]
1085 1083  
1086   - text_file = open(filename, "r")
1087   - header = text_file.readline().split('\t')
1088   - text_file.close()
  1084 + registration_coordinates = np.array(data[1:]).astype(np.float32)
  1085 + self.obj_fiducials = registration_coordinates[:, :3]
  1086 + self.obj_orients = registration_coordinates[:, 3:]
1089 1087  
1090   - self.obj_name = header[1]
1091   - self.obj_ref_mode = int(header[-1])
  1088 + self.obj_name = data[0][1]
  1089 + self.obj_ref_mode = int(data[0][-1])
1092 1090  
1093 1091 self.checktrack.Enable(1)
1094 1092 self.checktrack.SetValue(True)
... ... @@ -1098,7 +1096,7 @@ class ObjectRegistrationPanel(wx.Panel):
1098 1096 label=_("Object file successfully loaded"))
1099 1097 Publisher.sendMessage('Update track object state', flag=True, obj_name=self.obj_name)
1100 1098 Publisher.sendMessage('Change camera checkbox', status=False)
1101   - # wx.MessageBox(_("Object file successfully loaded"), _("Load"))
  1099 + wx.MessageBox(_("Object file successfully loaded"), _("InVesalius 3"))
1102 1100 except:
1103 1101 wx.MessageBox(_("Object registration file incompatible."), _("InVesalius 3"))
1104 1102 Publisher.sendMessage('Update status text in GUI', label="")
... ... @@ -1449,8 +1447,8 @@ class MarkersPanel(wx.Panel):
1449 1447 else:
1450 1448 self.nav_status = True
1451 1449  
1452   - def UpdateSeedCoordinates(self, root=None, affine_vtk=None, coord_offset=(0, 0, 0)):
1453   - self.current_seed = coord_offset
  1450 + def UpdateSeedCoordinates(self, root=None, affine_vtk=None, coord_offset=(0, 0, 0), coord_offset_w=(0, 0, 0)):
  1451 + self.current_seed = coord_offset_w
1454 1452  
1455 1453 def UpdateRobotCoordinates(self, coordinates_raw, markers_flag):
1456 1454 self.raw_target_robot = coordinates_raw[1], coordinates_raw[2]
... ... @@ -1735,7 +1733,7 @@ class TractographyPanel(wx.Panel):
1735 1733 default_colour = wx.SystemSettings_GetColour(wx.SYS_COLOUR_MENUBAR)
1736 1734 self.SetBackgroundColour(default_colour)
1737 1735  
1738   - self.affine = None
  1736 + self.affine = np.identity(4)
1739 1737 self.affine_vtk = None
1740 1738 self.trekker = None
1741 1739 self.n_tracts = const.N_TRACTS
... ... @@ -1992,7 +1990,6 @@ class TractographyPanel(wx.Panel):
1992 1990 self.nav_status = nav_status
1993 1991  
1994 1992 def OnLinkBrain(self, event=None):
1995   - Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
1996 1993 Publisher.sendMessage('Begin busy cursor')
1997 1994 mask_path = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, _("Import brain mask"))
1998 1995 img_path = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, _("Import T1 anatomical image"))
... ... @@ -2006,31 +2003,37 @@ class TractographyPanel(wx.Panel):
2006 2003 slic = sl.Slice()
2007 2004 prj_data = prj.Project()
2008 2005 matrix_shape = tuple(prj_data.matrix_shape)
  2006 + spacing = tuple(prj_data.spacing)
  2007 + img_shift = spacing[1] * (matrix_shape[1] - 1)
2009 2008 self.affine = slic.affine.copy()
2010   - self.affine[1, -1] -= matrix_shape[1]
  2009 + self.affine[1, -1] -= img_shift
2011 2010 self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine)
2012 2011  
2013   - try:
2014   - self.brain_peel = brain.Brain(img_path, mask_path, self.n_peels, self.affine_vtk)
2015   - self.brain_actor = self.brain_peel.get_actor(self.peel_depth, self.affine_vtk)
2016   - self.brain_actor.GetProperty().SetOpacity(self.brain_opacity)
2017   - Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor)
2018   - Publisher.sendMessage('Get peel centers and normals', centers=self.brain_peel.peel_centers,
2019   - normals=self.brain_peel.peel_normals)
2020   - Publisher.sendMessage('Get init locator', locator=self.brain_peel.locator)
2021   - self.checkpeeling.Enable(1)
2022   - self.checkpeeling.SetValue(True)
2023   - self.spin_opacity.Enable(1)
2024   - Publisher.sendMessage('Update status text in GUI', label=_("Brain model loaded"))
2025   - self.peel_loaded = True
2026   - Publisher.sendMessage('Update peel visualization', data= self.peel_loaded)
2027   - except:
2028   - wx.MessageBox(_("Unable to load brain mask."), _("InVesalius 3"))
  2012 + if mask_path and img_path:
  2013 + Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
  2014 + try:
  2015 + self.brain_peel = brain.Brain(img_path, mask_path, self.n_peels, self.affine_vtk)
  2016 + self.brain_actor = self.brain_peel.get_actor(self.peel_depth, self.affine_vtk)
  2017 + self.brain_actor.GetProperty().SetOpacity(self.brain_opacity)
  2018 +
  2019 + self.checkpeeling.Enable(1)
  2020 + self.checkpeeling.SetValue(True)
  2021 + self.spin_opacity.Enable(1)
  2022 + self.peel_loaded = True
  2023 +
  2024 + Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor)
  2025 + Publisher.sendMessage('Get peel centers and normals', centers=self.brain_peel.peel_centers,
  2026 + normals=self.brain_peel.peel_normals)
  2027 + Publisher.sendMessage('Get init locator', locator=self.brain_peel.locator)
  2028 + Publisher.sendMessage('Update status text in GUI', label=_("Brain model loaded"))
  2029 + Publisher.sendMessage('Update peel visualization', data= self.peel_loaded)
  2030 + except:
  2031 + Publisher.sendMessage('Update status text in GUI', label=_("Brain mask initialization failed."))
  2032 + wx.MessageBox(_("Unable to load brain mask."), _("InVesalius 3"))
2029 2033  
2030 2034 Publisher.sendMessage('End busy cursor')
2031 2035  
2032 2036 def OnLinkFOD(self, event=None):
2033   - Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
2034 2037 Publisher.sendMessage('Begin busy cursor')
2035 2038 filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import Trekker FOD"))
2036 2039 # Juuso
... ... @@ -2041,69 +2044,83 @@ class TractographyPanel(wx.Panel):
2041 2044 # FOD_path = 'Baran_FOD.nii'
2042 2045 # filename = os.path.join(data_dir, FOD_path)
2043 2046  
2044   - # if not self.affine_vtk:
2045   - # slic = sl.Slice()
2046   - # self.affine = slic.affine
2047   - # self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine)
2048   -
2049 2047 if not self.affine_vtk:
2050 2048 slic = sl.Slice()
2051 2049 prj_data = prj.Project()
2052 2050 matrix_shape = tuple(prj_data.matrix_shape)
  2051 + spacing = tuple(prj_data.spacing)
  2052 + img_shift = spacing[1] * (matrix_shape[1] - 1)
2053 2053 self.affine = slic.affine.copy()
2054   - self.affine[1, -1] -= matrix_shape[1]
  2054 + self.affine[1, -1] -= img_shift
2055 2055 self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine)
2056 2056  
2057   - # try:
2058   -
2059   - self.trekker = Trekker.initialize(filename.encode('utf-8'))
2060   - self.trekker, n_threads = dti.set_trekker_parameters(self.trekker, self.trekker_cfg)
2061   -
2062   - self.checktracts.Enable(1)
2063   - self.checktracts.SetValue(True)
2064   - self.view_tracts = True
2065   - Publisher.sendMessage('Update Trekker object', data=self.trekker)
2066   - Publisher.sendMessage('Update number of threads', data=n_threads)
2067   - Publisher.sendMessage('Update tracts visualization', data=1)
2068   - Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialized"))
2069   - # except:
2070   - # wx.MessageBox(_("Unable to initialize Trekker, check FOD and config files."), _("InVesalius 3"))
  2057 + if filename:
  2058 + Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
  2059 + try:
  2060 + self.trekker = Trekker.initialize(filename.encode('utf-8'))
  2061 + self.trekker, n_threads = dti.set_trekker_parameters(self.trekker, self.trekker_cfg)
  2062 +
  2063 + self.checktracts.Enable(1)
  2064 + self.checktracts.SetValue(True)
  2065 + self.view_tracts = True
  2066 +
  2067 + Publisher.sendMessage('Update Trekker object', data=self.trekker)
  2068 + Publisher.sendMessage('Update number of threads', data=n_threads)
  2069 + Publisher.sendMessage('Update tracts visualization', data=1)
  2070 + Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialized"))
  2071 + # except:
  2072 + # wx.MessageBox(_("Unable to initialize Trekker, check FOD and config files."), _("InVesalius 3"))
  2073 + except:
  2074 + Publisher.sendMessage('Update status text in GUI', label=_("Trekker initialization failed."))
  2075 + wx.MessageBox(_("Unable to load FOD."), _("InVesalius 3"))
2071 2076  
2072 2077 Publisher.sendMessage('End busy cursor')
2073 2078  
2074 2079 def OnLoadACT(self, event=None):
2075   - Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
2076   - Publisher.sendMessage('Begin busy cursor')
2077   - filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import anatomical labels"))
2078   - # Baran
2079   - # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\anat_reg_improve_20200609'
2080   - # act_path = 'Baran_trekkerACTlabels_inFODspace.nii'
2081   - # filename = os.path.join(data_dir, act_path)
2082   -
2083   - act_data = nb.squeeze_image(nb.load(filename))
2084   - act_data = nb.as_closest_canonical(act_data)
2085   - act_data.update_header()
2086   - act_data_arr = act_data.get_fdata()
  2080 + if self.trekker:
  2081 + Publisher.sendMessage('Begin busy cursor')
  2082 + filename = dlg.ShowImportOtherFilesDialog(const.ID_NIFTI_IMPORT, msg=_("Import anatomical labels"))
  2083 + # Baran
  2084 + # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\anat_reg_improve_20200609'
  2085 + # act_path = 'Baran_trekkerACTlabels_inFODspace.nii'
  2086 + # filename = os.path.join(data_dir, act_path)
  2087 +
  2088 + if not self.affine_vtk:
  2089 + slic = sl.Slice()
  2090 + prj_data = prj.Project()
  2091 + matrix_shape = tuple(prj_data.matrix_shape)
  2092 + spacing = tuple(prj_data.spacing)
  2093 + img_shift = spacing[1] * (matrix_shape[1] - 1)
  2094 + self.affine = slic.affine.copy()
  2095 + self.affine[1, -1] -= img_shift
  2096 + self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine)
2087 2097  
2088   - if not self.affine_vtk:
2089   - slic = sl.Slice()
2090   - prj_data = prj.Project()
2091   - matrix_shape = tuple(prj_data.matrix_shape)
2092   - self.affine = slic.affine.copy()
2093   - self.affine[1, -1] -= matrix_shape[1]
2094   - self.affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(self.affine)
2095   -
2096   - self.checkACT.Enable(1)
2097   - self.checkACT.SetValue(True)
2098   -
2099   - Publisher.sendMessage('Update ACT data', data=act_data_arr)
2100   - Publisher.sendMessage('Enable ACT', data=True)
2101   - # Publisher.sendMessage('Create grid', data=act_data_arr, affine=self.affine)
2102   - # Publisher.sendMessage('Update number of threads', data=n_threads)
2103   - # Publisher.sendMessage('Update tracts visualization', data=1)
2104   - Publisher.sendMessage('Update status text in GUI', label=_("Trekker ACT loaded"))
2105   -
2106   - Publisher.sendMessage('End busy cursor')
  2098 + try:
  2099 + Publisher.sendMessage('Update status text in GUI', label=_("Busy"))
  2100 + if filename:
  2101 + act_data = nb.squeeze_image(nb.load(filename))
  2102 + act_data = nb.as_closest_canonical(act_data)
  2103 + act_data.update_header()
  2104 + act_data_arr = act_data.get_fdata()
  2105 +
  2106 + self.checkACT.Enable(1)
  2107 + self.checkACT.SetValue(True)
  2108 +
  2109 + # ACT rules should be as follows:
  2110 + self.trekker.pathway_stop_at_entry(filename.encode('utf-8'), -1) # outside
  2111 + self.trekker.pathway_discard_if_ends_inside(filename.encode('utf-8'), 1) # wm
  2112 + self.trekker.pathway_discard_if_enters(filename.encode('utf-8'), 0) # csf
  2113 +
  2114 + Publisher.sendMessage('Update ACT data', data=act_data_arr)
  2115 + Publisher.sendMessage('Enable ACT', data=True)
  2116 + Publisher.sendMessage('Update status text in GUI', label=_("Trekker ACT loaded"))
  2117 + except:
  2118 + Publisher.sendMessage('Update status text in GUI', label=_("ACT initialization failed."))
  2119 + wx.MessageBox(_("Unable to load ACT."), _("InVesalius 3"))
  2120 +
  2121 + Publisher.sendMessage('End busy cursor')
  2122 + else:
  2123 + wx.MessageBox(_("Load FOD image before the ACT."), _("InVesalius 3"))
2107 2124  
2108 2125 def OnLoadParameters(self, event=None):
2109 2126 import json
... ... @@ -2146,10 +2163,13 @@ class TractographyPanel(wx.Panel):
2146 2163 # print("Running during navigation")
2147 2164 coord_flip = list(position[:3])
2148 2165 coord_flip[1] = -coord_flip[1]
2149   - dti.compute_tracts(self.trekker, coord_flip, self.affine, self.affine_vtk,
2150   - self.n_tracts)
  2166 + dti.compute_and_visualize_tracts(self.trekker, coord_flip, self.affine, self.affine_vtk,
  2167 + self.n_tracts)
2151 2168  
2152 2169 def OnCloseProject(self):
  2170 + self.trekker = None
  2171 + self.trekker_cfg = const.TREKKER_CONFIG
  2172 +
2153 2173 self.checktracts.SetValue(False)
2154 2174 self.checktracts.Enable(0)
2155 2175 self.checkpeeling.SetValue(False)
... ...
invesalius/navigation/navigation.py
  Diff comments   View file @46e9180
... ... @@ -27,7 +27,6 @@ import numpy as np
27 27 import invesalius.constants as const
28 28 import invesalius.project as prj
29 29 import invesalius.data.bases as db
30   -import invesalius.data.coordinates as dco
31 30 import invesalius.data.coregistration as dcr
32 31 import invesalius.data.serial_port_connection as spc
33 32 import invesalius.data.slice_ as sl
... ... @@ -98,12 +97,11 @@ class UpdateNavigationScene(threading.Thread):
98 97  
99 98 # use of CallAfter is mandatory otherwise crashes the wx interface
100 99 if self.view_tracts:
101   - bundle, affine_vtk, coord_offset = self.tracts_queue.get_nowait()
  100 + bundle, affine_vtk, coord_offset, coord_offset_w = self.tracts_queue.get_nowait()
102 101 #TODO: Check if possible to combine the Remove tracts with Update tracts in a single command
103 102 wx.CallAfter(Publisher.sendMessage, 'Remove tracts')
104   - wx.CallAfter(Publisher.sendMessage, 'Update tracts', root=bundle,
105   - affine_vtk=affine_vtk, coord_offset=coord_offset)
106   - # wx.CallAfter(Publisher.sendMessage, 'Update marker offset', coord_offset=coord_offset)
  103 + wx.CallAfter(Publisher.sendMessage, 'Update tracts', root=bundle, affine_vtk=affine_vtk,
  104 + coord_offset=coord_offset, coord_offset_w=coord_offset_w)
107 105 self.tracts_queue.task_done()
108 106  
109 107 if self.serial_port_enabled:
... ... @@ -118,10 +116,11 @@ class UpdateNavigationScene(threading.Thread):
118 116 wx.CallAfter(Publisher.sendMessage, 'Set cross focal point', position=coord)
119 117 wx.CallAfter(Publisher.sendMessage, 'Update raw coordinates', coordinates_raw=coordinates_raw, markers_flag=markers_flag)
120 118 wx.CallAfter(Publisher.sendMessage, 'Update slice viewer')
  119 + wx.CallAfter(Publisher.sendMessage, 'Sensor ID', markers_flag=markers_flag)
121 120  
122 121 if view_obj:
123 122 wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord)
124   - wx.CallAfter(Publisher.sendMessage, 'Update object arrow matrix',m_img=m_img, coord=coord, flag= self.peel_loaded)
  123 + wx.CallAfter(Publisher.sendMessage, 'Update object arrow matrix', m_img=m_img, coord=coord, flag= self.peel_loaded)
125 124 self.coord_queue.task_done()
126 125 # print('UpdateScene: done {}'.format(count))
127 126 # count += 1
... ... @@ -191,7 +190,7 @@ class Navigation():
191 190  
192 191 def UpdateSleep(self, sleep):
193 192 self.sleep_nav = sleep
194   - self.serial_port_connection.sleep_nav = sleep
  193 + # self.serial_port_connection.sleep_nav = sleep
195 194  
196 195 def UpdateSerialPort(self, serial_port_in_use, com_port=None, baud_rate=None):
197 196 self.serial_port_in_use = serial_port_in_use
... ... @@ -307,15 +306,22 @@ class Navigation():
307 306  
308 307 if self.view_tracts:
309 308 # initialize Trekker parameters
  309 + # TODO: This affine and affine_vtk are created 4 times. To improve, create a new affine object inside
  310 + # Slice() that contains the transformation with the img_shift. Rename it to avoid confusion to the
  311 + # affine, for instance it can be: affine_world_to_invesalius_vtk
310 312 slic = sl.Slice()
311 313 prj_data = prj.Project()
312 314 matrix_shape = tuple(prj_data.matrix_shape)
  315 + spacing = tuple(prj_data.spacing)
  316 + img_shift = spacing[1] * (matrix_shape[1] - 1)
313 317 affine = slic.affine.copy()
314   - affine[1, -1] -= matrix_shape[1]
  318 + affine[1, -1] -= img_shift
315 319 affine_vtk = vtk_utils.numpy_to_vtkMatrix4x4(affine)
  320 +
316 321 Publisher.sendMessage("Update marker offset state", create=True)
  322 +
317 323 self.trk_inp = self.trekker, affine, self.seed_offset, self.n_tracts, self.seed_radius,\
318   - self.n_threads, self.act_data, affine_vtk, matrix_shape[1]
  324 + self.n_threads, self.act_data, affine_vtk, img_shift
319 325 # print("Appending the tract computation thread!")
320 326 queues = [self.coord_tracts_queue, self.tracts_queue]
321 327 if self.enable_act:
... ...