From 0673940c840118edf2b74d57dd93d5dd648a61a0 Mon Sep 17 00:00:00 2001 From: Victor Hugo Souza Date: Thu, 20 Aug 2020 17:21:01 +0300 Subject: [PATCH] Refactor of thread management in neuronavigation (#242) --- invesalius/constants.py | 25 ++++++++++++++++++++++++- invesalius/control.py | 89 +++++++++++++++++++++++++++++------------------------------------------------------------ invesalius/data/bases.py | 244 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------ invesalius/data/coordinates.py | 138 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------- invesalius/data/coregistration.py | 911 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/imagedata_utils.py | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius/data/record_coords.py | 8 ++++++-- invesalius/data/slice_.py | 27 +++++++++++++++++++++++++++ invesalius/data/styles.py | 124 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---- invesalius/data/surface.py | 3 +++ invesalius/data/trackers.py | 3 +++ invesalius/data/tractography.py | 471 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius/data/trigger.py | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- invesalius/data/viewer_slice.py | 18 +++++++++--------- invesalius/data/viewer_volume.py | 310 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------ invesalius/data/vtk_utils.py | 18 ++++++++++++++++++ invesalius/gui/dialogs.py | 265 ++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/gui/frame.py | 7 ++++--- invesalius/gui/task_navigator.py | 902 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/project.py | 3 ++- 20 files changed, 2711 insertions(+), 1011 deletions(-) create mode 100644 invesalius/data/tractography.py diff --git a/invesalius/constants.py b/invesalius/constants.py index 95b814c..53e1ac2 100644 --- a/invesalius/constants.py +++ b/invesalius/constants.py @@ -19,6 +19,7 @@ import os.path import platform +import psutil import sys import wx import itertools @@ -521,6 +522,11 @@ ID_GOTO_COORD = wx.NewId() ID_MANUAL_WWWL = wx.NewId() +# Tractography with Trekker +ID_TREKKER_MASK = wx.NewId() +ID_TREKKER_IMG = wx.NewId() +ID_TREKKER_FOD = wx.NewId() + #--------------------------------------------------------- STATE_DEFAULT = 1000 STATE_WL = 1001 @@ -545,6 +551,7 @@ SLICE_STATE_REMOVE_MASK_PARTS = 3012 SLICE_STATE_SELECT_MASK_PARTS = 3013 SLICE_STATE_FFILL_SEGMENTATION = 3014 SLICE_STATE_CROP_MASK = 3015 +SLICE_STATE_TRACTS = 3016 VOLUME_STATE_SEED = 2001 # STATE_LINEAR_MEASURE = 3001 @@ -559,7 +566,7 @@ TOOL_STATES = [STATE_WL, STATE_SPIN, STATE_ZOOM, TOOL_SLICE_STATES = [SLICE_STATE_CROSS, SLICE_STATE_SCROLL, - SLICE_STATE_REORIENT] + SLICE_STATE_REORIENT, SLICE_STATE_TRACTS] SLICE_STYLES = TOOL_STATES + TOOL_SLICE_STATES @@ -651,6 +658,8 @@ BOOLEAN_XOR = 4 #------------ Navigation defaults ------------------- +MARKER_COLOUR = (1.0, 1.0, 0.) +MARKER_SIZE = 2 SELECT = 0 MTC = 1 FASTRAK = 2 @@ -738,3 +747,17 @@ COIL_COORD_THRESHOLD = 3 TIMESTAMP = 2.0 CAM_MODE = True + +# Tractography visualization +N_TRACTS = 100 +PEEL_DEPTH = 5 +MAX_PEEL_DEPTH = 30 +SEED_OFFSET = 15 +SEED_RADIUS = 1.5 +SLEEP_NAVIGATION = 0.1 +BRAIN_OPACITY = 0.5 +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} diff --git a/invesalius/control.py b/invesalius/control.py index f288077..e2d602e 100644 --- a/invesalius/control.py +++ b/invesalius/control.py @@ -32,6 +32,7 @@ import invesalius.data.mask as msk import invesalius.data.measures as measures import invesalius.data.slice_ as sl import invesalius.data.surface as srf +import invesalius.data.transformations as tr import invesalius.data.volume as volume import invesalius.gui.dialogs as dialog import invesalius.project as prj @@ -57,7 +58,6 @@ class Controller(): def __init__(self, frame): self.surface_manager = srf.SurfaceManager() self.volume = volume.Volume() - self.plugin_manager = plugins.PluginManager() self.__bind_events() self.frame = frame self.progress_dialog = None @@ -76,12 +76,9 @@ class Controller(): Publisher.sendMessage('Load Preferences') - self.plugin_manager.find_plugins() - def __bind_events(self): Publisher.subscribe(self.OnImportMedicalImages, 'Import directory') Publisher.subscribe(self.OnImportGroup, 'Import group') - Publisher.subscribe(self.OnImportFolder, 'Import folder') Publisher.subscribe(self.OnShowDialogImportDirectory, 'Show import directory dialog') Publisher.subscribe(self.OnShowDialogImportOtherFiles, @@ -123,8 +120,6 @@ class Controller(): Publisher.subscribe(self.Send_affine, 'Get affine matrix') - Publisher.subscribe(self.create_project_from_matrix, 'Create project from matrix') - def SetBitmapSpacing(self, spacing): proj = prj.Project() proj.spacing = spacing @@ -336,6 +331,10 @@ class Controller(): self.Slice.window_level = proj.level self.Slice.window_width = proj.window + if proj.affine: + self.Slice.affine = np.asarray(proj.affine).reshape(4, 4) + else: + self.Slice.affine = None Publisher.sendMessage('Update threshold limits list', threshold_range=proj.threshold_range) @@ -504,32 +503,7 @@ class Controller(): self.LoadProject() Publisher.sendMessage("Enable state project", state=True) - def OnImportFolder(self, folder): - Publisher.sendMessage('Begin busy cursor') - folder = os.path.abspath(folder) - - proj = prj.Project() - proj.load_from_folder(folder) - - self.Slice = sl.Slice() - self.Slice._open_image_matrix(proj.matrix_filename, - tuple(proj.matrix_shape), - proj.matrix_dtype) - - self.Slice.window_level = proj.level - self.Slice.window_width = proj.window - - Publisher.sendMessage('Update threshold limits list', - threshold_range=proj.threshold_range) - - session = ses.Session() - filename = proj.name+".inv3" - filename = filename.replace("/", "") #Fix problem case other/Skull_DICOM - dirpath = session.CreateProject(filename) - self.LoadProject() - Publisher.sendMessage("Enable state project", state=True) - - Publisher.sendMessage('End busy cursor') + #------------------------------------------------------------------------------------- def LoadProject(self): proj = prj.Project() @@ -688,6 +662,9 @@ class Controller(): proj.matrix_dtype = matrix.dtype.name proj.matrix_filename = matrix_filename + if self.affine is not None: + proj.affine = self.affine.tolist() + # Orientation must be CORONAL in order to as_closes_canonical and # swap axis in img2memmap to work in a standardized way. # TODO: Create standard import image for all acquisition orientations @@ -700,7 +677,6 @@ class Controller(): proj.level = self.Slice.window_level proj.threshold_range = int(matrix.min()), int(matrix.max()) proj.spacing = self.Slice.spacing - proj.affine = self.affine.tolist() ###### session = ses.Session() @@ -872,11 +848,13 @@ class Controller(): def OnOpenOtherFiles(self, filepath): filepath = utils.decode(filepath, const.FS_ENCODE) if not(filepath) == None: - name = os.path.basename(filepath).split(".")[0] + name = filepath.rpartition('\\')[-1].split('.') + group = oth.ReadOthers(filepath) + if group: matrix, matrix_filename = self.OpenOtherFiles(group) - self.CreateOtherProject(name, matrix, matrix_filename) + self.CreateOtherProject(str(name[0]), matrix, matrix_filename) self.LoadProject() Publisher.sendMessage("Enable state project", state=True) else: @@ -981,10 +959,10 @@ class Controller(): self.matrix, scalar_range, self.filename = image_utils.img2memmap(group) hdr = group.header - if group.affine.any(): - self.affine = group.affine - Publisher.sendMessage('Update affine matrix', - affine=self.affine, status=True) + # if group.affine.any(): + # self.affine = group.affine + # Publisher.sendMessage('Update affine matrix', + # affine=self.affine, status=True) hdr.set_data_dtype('int16') dims = hdr.get_zooms() dimsf = tuple([float(s) for s in dims]) @@ -1000,6 +978,18 @@ class Controller(): self.Slice.window_level = wl self.Slice.window_width = ww + if group.affine.any(): + # TODO: replace the inverse of the affine by the actual affine in the whole code + # remove scaling factor for non-unitary voxel dimensions + # self.affine = image_utils.world2invspace(affine=group.affine) + scale, shear, angs, trans, persp = tr.decompose_matrix(group.affine) + self.affine = np.linalg.inv(tr.compose_matrix(scale=None, shear=shear, + angles=angs, translate=trans, perspective=persp)) + # print("repos_img: {}".format(repos_img)) + self.Slice.affine = self.affine + Publisher.sendMessage('Update affine matrix', + affine=self.affine, status=True) + scalar_range = int(scalar_range[0]), int(scalar_range[1]) Publisher.sendMessage('Update threshold limits list', threshold_range=scalar_range) @@ -1066,24 +1056,3 @@ class Controller(): def ApplyReorientation(self): self.Slice.apply_reorientation() - - def start_new_inv_instance(self, image, name, spacing, modality, orientation, window_width, window_level): - p = prj.Project() - project_folder = tempfile.mkdtemp() - p.create_project_file(name, spacing, modality, orientation, window_width, window_level, image, folder=project_folder) - err_msg = '' - try: - sp = subprocess.Popen([sys.executable, sys.argv[0], '--import-folder', project_folder], - stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=os.getcwd()) - except Exception as err: - err_msg = str(err) - else: - try: - if sp.wait(2): - err_msg = sp.stderr.read().decode('utf8') - sp.terminate() - except subprocess.TimeoutExpired: - pass - - if err_msg: - dialog.MessageBox(None, "It was not possible to launch new instance of InVesalius3 dsfa dfdsfa sdfas fdsaf asdfasf dsaa", err_msg) diff --git a/invesalius/data/bases.py b/invesalius/data/bases.py index e43a0f5..de32b67 100644 --- a/invesalius/data/bases.py +++ b/invesalius/data/bases.py @@ -1,4 +1,3 @@ -from math import sqrt, pi import numpy as np import invesalius.data.coordinates as dco import invesalius.data.transformations as tr @@ -23,8 +22,7 @@ def angle_calculation(ap_axis, coil_axis): def base_creation_old(fiducials): """ - Calculate the origin and matrix for coordinate system - transformation. + Calculate the origin and matrix for coordinate system transformation. q: origin of coordinate system g1, g2, g3: orthogonal vectors of coordinate system @@ -49,9 +47,9 @@ def base_creation_old(fiducials): g3 = np.cross(g2, g1) - g1 = g1/sqrt(np.dot(g1, g1)) - g2 = g2/sqrt(np.dot(g2, g2)) - g3 = g3/sqrt(np.dot(g3, g3)) + g1 = g1/np.sqrt(np.dot(g1, g1)) + g2 = g2/np.sqrt(np.dot(g2, g2)) + g3 = g3/np.sqrt(np.dot(g3, g3)) m = np.matrix([[g1[0], g1[1], g1[2]], [g2[0], g2[1], g2[2]], @@ -79,7 +77,7 @@ def base_creation(fiducials): sub1 = p2 - p1 sub2 = p3 - p1 - lamb = (sub1[0]*sub2[0]+sub1[1]*sub2[1]+sub1[2]*sub2[2])/np.dot(sub1, sub1) + lamb = np.dot(sub1, sub2)/np.dot(sub1, sub1) q = p1 + lamb*sub1 g1 = p3 - q @@ -90,20 +88,52 @@ def base_creation(fiducials): g3 = np.cross(g1, g2) - g1 = g1/sqrt(np.dot(g1, g1)) - g2 = g2/sqrt(np.dot(g2, g2)) - g3 = g3/sqrt(np.dot(g3, g3)) - - m = np.matrix([[g1[0], g2[0], g3[0]], - [g1[1], g2[1], g3[1]], - [g1[2], g2[2], g3[2]]]) - - m_inv = m.I - - return m, q, m_inv - - -def calculate_fre(fiducials, minv, n, q, o): + g1 = g1/np.sqrt(np.dot(g1, g1)) + g2 = g2/np.sqrt(np.dot(g2, g2)) + g3 = g3/np.sqrt(np.dot(g3, g3)) + + m = np.zeros([3, 3]) + m[:, 0] = g1/np.sqrt(np.dot(g1, g1)) + m[:, 1] = g2/np.sqrt(np.dot(g2, g2)) + m[:, 2] = g3/np.sqrt(np.dot(g3, g3)) + + return m, q + + +# def calculate_fre(fiducials, minv, n, q, o): +# """ +# Calculate the Fiducial Registration Error for neuronavigation. +# +# :param fiducials: array of 6 rows (image and tracker fiducials) and 3 columns (x, y, z) with coordinates +# :param minv: inverse matrix given by base creation +# :param n: base change matrix given by base creation +# :param q: origin of first base +# :param o: origin of second base +# :return: float number of fiducial registration error +# """ +# +# img = np.zeros([3, 3]) +# dist = np.zeros([3, 1]) +# +# q1 = np.mat(q).reshape(3, 1) +# q2 = np.mat(o).reshape(3, 1) +# +# p1 = np.mat(fiducials[3, :]).reshape(3, 1) +# p2 = np.mat(fiducials[4, :]).reshape(3, 1) +# p3 = np.mat(fiducials[5, :]).reshape(3, 1) +# +# img[0, :] = np.asarray((q1 + (minv * n) * (p1 - q2)).reshape(1, 3)) +# img[1, :] = np.asarray((q1 + (minv * n) * (p2 - q2)).reshape(1, 3)) +# img[2, :] = np.asarray((q1 + (minv * n) * (p3 - q2)).reshape(1, 3)) +# +# dist[0] = np.sqrt(np.sum(np.power((img[0, :] - fiducials[0, :]), 2))) +# dist[1] = np.sqrt(np.sum(np.power((img[1, :] - fiducials[1, :]), 2))) +# dist[2] = np.sqrt(np.sum(np.power((img[2, :] - fiducials[2, :]), 2))) +# +# return float(np.sqrt(np.sum(dist ** 2) / 3)) + + +def calculate_fre_m(fiducials): """ Calculate the Fiducial Registration Error for neuronavigation. @@ -115,19 +145,29 @@ def calculate_fre(fiducials, minv, n, q, o): :return: float number of fiducial registration error """ + m, q1, minv = base_creation_old(fiducials[:3, :]) + n, q2, ninv = base_creation_old(fiducials[3:, :]) + + # TODO: replace the old by the new base creation + # the values differ greatly if FRE is computed using the old or new base_creation + # check the reason for the difference, because they should be the same + # m, q1 = base_creation(fiducials[:3, :]) + # n, q2 = base_creation(fiducials[3:, :]) + # minv = np.linalg.inv(m) + img = np.zeros([3, 3]) dist = np.zeros([3, 1]) - q1 = np.mat(q).reshape(3, 1) - q2 = np.mat(o).reshape(3, 1) + q1 = q1.reshape(3, 1) + q2 = q2.reshape(3, 1) - p1 = np.mat(fiducials[3, :]).reshape(3, 1) - p2 = np.mat(fiducials[4, :]).reshape(3, 1) - p3 = np.mat(fiducials[5, :]).reshape(3, 1) + p1 = fiducials[3, :].reshape(3, 1) + p2 = fiducials[4, :].reshape(3, 1) + p3 = fiducials[5, :].reshape(3, 1) - img[0, :] = np.asarray((q1 + (minv * n) * (p1 - q2)).reshape(1, 3)) - img[1, :] = np.asarray((q1 + (minv * n) * (p2 - q2)).reshape(1, 3)) - img[2, :] = np.asarray((q1 + (minv * n) * (p3 - q2)).reshape(1, 3)) + img[0, :] = (q1 + (minv @ n) * (p1 - q2)).reshape(1, 3) + img[1, :] = (q1 + (minv @ n) * (p2 - q2)).reshape(1, 3) + img[2, :] = (q1 + (minv @ n) * (p3 - q2)).reshape(1, 3) dist[0] = np.sqrt(np.sum(np.power((img[0, :] - fiducials[0, :]), 2))) dist[1] = np.sqrt(np.sum(np.power((img[1, :] - fiducials[1, :]), 2))) @@ -136,44 +176,68 @@ def calculate_fre(fiducials, minv, n, q, o): return float(np.sqrt(np.sum(dist ** 2) / 3)) -def flip_x(point): - """ - Flip coordinates of a vector according to X axis - Coronal Images do not require this transformation - 1 tested - and for this case, at navigation, the z axis is inverted - - It's necessary to multiply the z coordinate by (-1). Possibly - because the origin of coordinate system of imagedata is - located in superior left corner and the origin of VTK scene coordinate - system (polygonal surface) is in the interior left corner. Second - possibility is the order of slice stacking - - :param point: list of coordinates x, y and z - :return: flipped coordinates - """ - - # TODO: check if the Flip function is related to the X or Y axis - - point = np.matrix(point + (0,)) - point[0, 2] = -point[0, 2] - - m_rot = np.matrix([[1.0, 0.0, 0.0, 0.0], - [0.0, -1.0, 0.0, 0.0], - [0.0, 0.0, -1.0, 0.0], - [0.0, 0.0, 0.0, 1.0]]) - m_trans = np.matrix([[1.0, 0, 0, -point[0, 0]], - [0.0, 1.0, 0, -point[0, 1]], - [0.0, 0.0, 1.0, -point[0, 2]], - [0.0, 0.0, 0.0, 1.0]]) - m_trans_return = np.matrix([[1.0, 0, 0, point[0, 0]], - [0.0, 1.0, 0, point[0, 1]], - [0.0, 0.0, 1.0, point[0, 2]], - [0.0, 0.0, 0.0, 1.0]]) - - point_rot = point*m_trans*m_rot*m_trans_return - x, y, z = point_rot.tolist()[0][:3] - - return x, y, z +# def flip_x(point): +# """ +# Flip coordinates of a vector according to X axis +# Coronal Images do not require this transformation - 1 tested +# and for this case, at navigation, the z axis is inverted +# +# It's necessary to multiply the z coordinate by (-1). Possibly +# because the origin of coordinate system of imagedata is +# located in superior left corner and the origin of VTK scene coordinate +# system (polygonal surface) is in the interior left corner. Second +# possibility is the order of slice stacking +# +# :param point: list of coordinates x, y and z +# :return: flipped coordinates +# """ +# +# # TODO: check if the Flip function is related to the X or Y axis +# +# point = np.matrix(point + (0,)) +# point[0, 2] = -point[0, 2] +# +# m_rot = np.matrix([[1.0, 0.0, 0.0, 0.0], +# [0.0, -1.0, 0.0, 0.0], +# [0.0, 0.0, -1.0, 0.0], +# [0.0, 0.0, 0.0, 1.0]]) +# m_trans = np.matrix([[1.0, 0, 0, -point[0, 0]], +# [0.0, 1.0, 0, -point[0, 1]], +# [0.0, 0.0, 1.0, -point[0, 2]], +# [0.0, 0.0, 0.0, 1.0]]) +# m_trans_return = np.matrix([[1.0, 0, 0, point[0, 0]], +# [0.0, 1.0, 0, point[0, 1]], +# [0.0, 0.0, 1.0, point[0, 2]], +# [0.0, 0.0, 0.0, 1.0]]) +# +# point_rot = point*m_trans*m_rot*m_trans_return +# x, y, z = point_rot.tolist()[0][:3] +# +# return x, y, z + + +# def flip_x_m(point): +# """ +# Rotate coordinates of a vector by pi around X axis in static reference frame. +# +# InVesalius also require to multiply the z coordinate by (-1). Possibly +# because the origin of coordinate system of imagedata is +# located in superior left corner and the origin of VTK scene coordinate +# system (polygonal surface) is in the interior left corner. Second +# possibility is the order of slice stacking +# +# :param point: list of coordinates x, y and z +# :return: rotated coordinates +# """ +# +# point_4 = np.hstack((point, 1.)).reshape([4, 1]) +# point_4[2, 0] = -point_4[2, 0] +# +# m_rot = tr.euler_matrix(pi, 0, 0) +# +# point_rot = m_rot @ point_4 +# +# return point_rot[0, 0], point_rot[1, 0], point_rot[2, 0] def flip_x_m(point): @@ -190,14 +254,14 @@ def flip_x_m(point): :return: rotated coordinates """ - point_4 = np.hstack((point, 1.)).reshape([4, 1]) + point_4 = np.hstack((point, 1.)).reshape(4, 1) point_4[2, 0] = -point_4[2, 0] - m_rot = np.asmatrix(tr.euler_matrix(pi, 0, 0)) + m_rot = tr.euler_matrix(np.pi, 0, 0) - point_rot = m_rot*point_4 + point_rot = m_rot @ point_4 - return point_rot[0, 0], point_rot[1, 0], point_rot[2, 0] + return point_rot def object_registration(fiducials, orients, coord_raw, m_change): @@ -224,48 +288,48 @@ def object_registration(fiducials, orients, coord_raw, m_change): fids_raw[ic, :] = dco.dynamic_reference_m2(coords[ic, :], coords[3, :])[:3] # compute initial alignment of probe fixed in the object in source frame - t_s0_raw = np.asmatrix(tr.translation_matrix(coords[3, :3])) - r_s0_raw = np.asmatrix(tr.euler_matrix(np.radians(coords[3, 3]), np.radians(coords[3, 4]), - np.radians(coords[3, 5]), 'rzyx')) - s0_raw = np.asmatrix(tr.concatenate_matrices(t_s0_raw, r_s0_raw)) + t_s0_raw = tr.translation_matrix(coords[3, :3]) + r_s0_raw = tr.euler_matrix(np.radians(coords[3, 3]), np.radians(coords[3, 4]), + np.radians(coords[3, 5]), 'rzyx') + s0_raw = tr.concatenate_matrices(t_s0_raw, r_s0_raw) # compute change of basis for object fiducials in source frame - base_obj_raw, q_obj_raw, base_inv_obj_raw = base_creation(fids_raw[:3, :3]) - r_obj_raw = np.asmatrix(np.identity(4)) + base_obj_raw, q_obj_raw = base_creation(fids_raw[:3, :3]) + r_obj_raw = np.identity(4) r_obj_raw[:3, :3] = base_obj_raw[:3, :3] - t_obj_raw = np.asmatrix(tr.translation_matrix(q_obj_raw)) - m_obj_raw = np.asmatrix(tr.concatenate_matrices(t_obj_raw, r_obj_raw)) + t_obj_raw = tr.translation_matrix(q_obj_raw) + m_obj_raw = tr.concatenate_matrices(t_obj_raw, r_obj_raw) for ic in range(0, 4): if coord_raw.any(): # compute object fiducials in reference frame fids_dyn[ic, :] = dco.dynamic_reference_m2(coords[ic, :], coord_raw[1, :]) - fids_dyn[ic, 2] = -fids_dyn[ic, 2] else: # compute object fiducials in source frame fids_dyn[ic, :] = coords[ic, :] + fids_dyn[ic, 2] = -fids_dyn[ic, 2] # compute object fiducials in vtk head frame a, b, g = np.radians(fids_dyn[ic, 3:]) T_p = tr.translation_matrix(fids_dyn[ic, :3]) R_p = tr.euler_matrix(a, b, g, 'rzyx') - M_p = np.asmatrix(tr.concatenate_matrices(T_p, R_p)) - M_img = np.asmatrix(m_change) * M_p + M_p = tr.concatenate_matrices(T_p, R_p) + M_img = m_change @ M_p angles_img = np.degrees(np.asarray(tr.euler_from_matrix(M_img, 'rzyx'))) - coord_img = np.asarray(flip_x_m(tr.translation_from_matrix(M_img))) + coord_img = flip_x_m(tr.translation_from_matrix(M_img)) - fids_img[ic, :] = np.hstack((coord_img, angles_img)) + fids_img[ic, :] = np.hstack((coord_img[:3, 0], angles_img)) # compute object base change in vtk head frame - base_obj_img, q_obj_img, base_inv_obj_img = base_creation(fids_img[:3, :3]) - r_obj_img = np.asmatrix(np.identity(4)) + base_obj_img, _ = base_creation(fids_img[:3, :3]) + r_obj_img = np.identity(4) r_obj_img[:3, :3] = base_obj_img[:3, :3] # compute initial alignment of probe fixed in the object in reference (or static) frame - s0_trans_dyn = np.asmatrix(tr.translation_matrix(fids_dyn[3, :3])) - s0_rot_dyn = np.asmatrix(tr.euler_matrix(np.radians(fids_dyn[3, 3]), np.radians(fids_dyn[3, 4]), - np.radians(fids_dyn[3, 5]), 'rzyx')) - s0_dyn = np.asmatrix(tr.concatenate_matrices(s0_trans_dyn, s0_rot_dyn)) + s0_trans_dyn = tr.translation_matrix(fids_dyn[3, :3]) + s0_rot_dyn = tr.euler_matrix(np.radians(fids_dyn[3, 3]), np.radians(fids_dyn[3, 4]), + np.radians(fids_dyn[3, 5]), 'rzyx') + s0_dyn = tr.concatenate_matrices(s0_trans_dyn, s0_rot_dyn) return t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img diff --git a/invesalius/data/coordinates.py b/invesalius/data/coordinates.py index 201ed2d..0bedb79 100644 --- a/invesalius/data/coordinates.py +++ b/invesalius/data/coordinates.py @@ -20,6 +20,7 @@ from math import sin, cos import numpy as np +import invesalius.data.bases as db import invesalius.data.transformations as tr import invesalius.constants as const @@ -74,7 +75,7 @@ def PolarisCoord(trck_init, trck_id, ref_mode): coord3 = np.hstack((trans_obj, angles_obj)) coord = np.vstack([coord1, coord2, coord3]) - Publisher.sendMessage('Sensors ID', probe_id=trck.probeID, ref_id=trck.refID, obj_id=trck.objID) + # Publisher.sendMessage('Sensors ID', probe_id=trck.probeID, ref_id=trck.refID, obj_id=trck.objID) return coord @@ -232,19 +233,45 @@ def DebugCoord(trk_init, trck_id, ref_mode): :return: six coordinates x, y, z, alfa, beta and gama """ - sleep(0.05) + # Started to take a more reasonable, limited random coordinate generator based on + # the collected fiducials, but it is more complicated than this. It should account for the + # dynamic reference computation + # if trk_init: + # fiducials = trk_init[3:, :] + # fids_max = fiducials.max(axis=0) + # fids_min = fiducials.min(axis=0) + # fids_lim = np.hstack((fids_min[np.newaxis, :].T, fids_max[np.newaxis, :].T)) + # + # dx = fids_max[] + # dt = [-180, 180] + # + # else: - coord1 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), - uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + dx = [-70, 70] + dt = [-180, 180] - coord2 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), - uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + coord1 = np.array([uniform(*dx), uniform(*dx), uniform(*dx), + uniform(*dt), uniform(*dt), uniform(*dt)]) + coord2 = np.array([uniform(*dx), uniform(*dx), uniform(*dx), + uniform(*dt), uniform(*dt), uniform(*dt)]) + coord3 = np.array([uniform(*dx), uniform(*dx), uniform(*dx), + uniform(*dt), uniform(*dt), uniform(*dt)]) + coord4 = np.array([uniform(*dx), uniform(*dx), uniform(*dx), + uniform(*dt), uniform(*dt), uniform(*dt)]) - coord3 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), - uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + sleep(0.05) - coord4 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), - uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + # coord1 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), + # uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + # + # coord2 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), + # uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + # + # coord3 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), + # uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) + # + # coord4 = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200), + # uniform(-180.0, 180.0), uniform(-180.0, 180.0), uniform(-180.0, 180.0)]) Publisher.sendMessage('Sensors ID', probe_id=int(uniform(0, 5)), ref_id=int(uniform(0, 5)), obj_id=int(uniform(0, 5))) @@ -297,19 +324,20 @@ def dynamic_reference_m(probe, reference): :param reference: sensor two defined as reference :return: rotated and translated coordinates """ - a, b, g = np.radians(reference[3:6]) - T = tr.translation_matrix(reference[:3]) - R = tr.euler_matrix(a, b, g, 'rzyx') - M = np.asmatrix(tr.concatenate_matrices(T, R)) - # M = tr.compose_matrix(angles=np.radians(reference[3:6]), translate=reference[:3]) - # print M - probe_4 = np.vstack((np.asmatrix(probe[:3]).reshape([3, 1]), 1.)) - coord_rot = M.I * probe_4 - coord_rot = np.squeeze(np.asarray(coord_rot)) + trans = tr.translation_matrix(reference[:3]) + rot = tr.euler_matrix(a, b, g, 'rzyx') + affine = tr.concatenate_matrices(trans, rot) + probe_4 = np.vstack((probe[:3].reshape([3, 1]), 1.)) + coord_rot = np.linalg.inv(affine) @ probe_4 + # minus sign to the z coordinate + coord_rot[2, 0] = -coord_rot[2, 0] + coord_rot = coord_rot[:3, 0].tolist() + coord_rot.extend(probe[3:]) + + return coord_rot - return coord_rot[0], coord_rot[1], -coord_rot[2], probe[3], probe[4], probe[5] def dynamic_reference_m2(probe, reference): """ @@ -331,70 +359,18 @@ def dynamic_reference_m2(probe, reference): T_p = tr.translation_matrix(probe[:3]) R = tr.euler_matrix(a, b, g, 'rzyx') R_p = tr.euler_matrix(a_p, b_p, g_p, 'rzyx') - M = np.asmatrix(tr.concatenate_matrices(T, R)) - M_p = np.asmatrix(tr.concatenate_matrices(T_p, R_p)) - # M = tr.compose_matrix(angles=np.radians(reference[3:6]), translate=reference[:3]) - # print M + M = tr.concatenate_matrices(T, R) + M_p = tr.concatenate_matrices(T_p, R_p) - M_dyn = M.I * M_p + M_dyn = np.linalg.inv(M) @ M_p al, be, ga = tr.euler_from_matrix(M_dyn, 'rzyx') coord_rot = tr.translation_from_matrix(M_dyn) coord_rot = np.squeeze(coord_rot) - # probe_4 = np.vstack((np.asmatrix(probe[:3]).reshape([3, 1]), 1.)) - # coord_rot_test = M.I * probe_4 - # coord_rot_test = np.squeeze(np.asarray(coord_rot_test)) - # - # print "coord_rot: ", coord_rot - # print "coord_rot_test: ", coord_rot_test - # print "test: ", np.allclose(coord_rot, coord_rot_test[:3]) - return coord_rot[0], coord_rot[1], coord_rot[2], np.degrees(al), np.degrees(be), np.degrees(ga) -# def dynamic_reference_m3(probe, reference): -# """ -# Apply dynamic reference correction to probe coordinates. Uses the alpha, beta and gama -# rotation angles of reference to rotate the probe coordinate and returns the x, y, z -# difference between probe and reference. Angles sequences and equation was extracted from -# Polhemus manual and Attitude matrix in Wikipedia. -# General equation is: -# coord = Mrot * (probe - reference) -# :param probe: sensor one defined as probe -# :param reference: sensor two defined as reference -# :return: rotated and translated coordinates -# """ -# -# a, b, g = np.radians(reference[3:6]) -# a_p, b_p, g_p = np.radians(probe[3:6]) -# -# T = tr.translation_matrix(reference[:3]) -# T_p = tr.translation_matrix(probe[:3]) -# R = tr.euler_matrix(a, b, g, 'rzyx') -# R_p = tr.euler_matrix(a_p, b_p, g_p, 'rzyx') -# M = np.asmatrix(tr.concatenate_matrices(T, R)) -# M_p = np.asmatrix(tr.concatenate_matrices(T_p, R_p)) -# # M = tr.compose_matrix(angles=np.radians(reference[3:6]), translate=reference[:3]) -# # print M -# -# M_dyn = M.I * M_p -# -# # al, be, ga = tr.euler_from_matrix(M_dyn, 'rzyx') -# # coord_rot = tr.translation_from_matrix(M_dyn) -# # -# # coord_rot = np.squeeze(coord_rot) -# -# # probe_4 = np.vstack((np.asmatrix(probe[:3]).reshape([3, 1]), 1.)) -# # coord_rot_test = M.I * probe_4 -# # coord_rot_test = np.squeeze(np.asarray(coord_rot_test)) -# # -# # print "coord_rot: ", coord_rot -# # print "coord_rot_test: ", coord_rot_test -# # print "test: ", np.allclose(coord_rot, coord_rot_test[:3]) -# -# return M_dyn - def str2float(data): """ @@ -414,3 +390,15 @@ def str2float(data): data = [float(s) for s in data[1:len(data)]] return data + + +def offset_coordinate(p_old, norm_vec, offset): + """ + Translate the coordinates of a point along a vector + :param p_old: (x, y, z) array with current point coordinates + :param norm_vec: (vx, vy, vz) array with normal vector coordinates + :param offset: double representing the magnitude of offset + :return: (x_new, y_new, z_new) array of offset coordinates + """ + p_offset = p_old - offset * norm_vec + return p_offset diff --git a/invesalius/data/coregistration.py b/invesalius/data/coregistration.py index 4f9e18a..df42805 100644 --- a/invesalius/data/coregistration.py +++ b/invesalius/data/coregistration.py @@ -17,347 +17,620 @@ # detalhes. #-------------------------------------------------------------------------- +import numpy as np +import queue import threading from time import sleep -from numpy import asmatrix, mat, degrees, radians, identity -import wx -from pubsub import pub as Publisher - import invesalius.data.coordinates as dco import invesalius.data.transformations as tr -# TODO: Optimize navigation thread. Remove the infinite loop and optimize sleep. - -class CoregistrationStatic(threading.Thread): - """ - Thread to update the coordinates with the fiducial points - co-registration method while the Navigation Button is pressed. - Sleep function in run method is used to avoid blocking GUI and - for better real-time navigation +# TODO: Replace the use of degrees by radians in every part of the navigation pipeline + +def object_marker_to_center(coord_raw, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw): + """Translate and rotate the raw coordinate given by the tracking device to the reference system created during + the object registration. + + :param coord_raw: Coordinates returned by the tracking device + :type coord_raw: numpy.ndarray + :param obj_ref_mode: + :type obj_ref_mode: int + :param t_obj_raw: + :type t_obj_raw: numpy.ndarray + :param s0_raw: + :type s0_raw: numpy.ndarray + :param r_s0_raw: rotation transformation from marker to object basis + :type r_s0_raw: numpy.ndarray + :return: 4 x 4 numpy double array + :rtype: numpy.ndarray """ - def __init__(self, coreg_data, nav_id, trck_info): - threading.Thread.__init__(self) - self.coreg_data = coreg_data - self.nav_id = nav_id - self.trck_info = trck_info - self._pause_ = False - self.start() - - def stop(self): - self._pause_ = True - - def run(self): - # m_change = self.coreg_data[0] - # obj_ref_mode = self.coreg_data[2] - # - # trck_init = self.trck_info[0] - # trck_id = self.trck_info[1] - # trck_mode = self.trck_info[2] - - m_change, obj_ref_mode = self.coreg_data - trck_init, trck_id, trck_mode = self.trck_info - - while self.nav_id: - coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) - - psi, theta, phi = coord_raw[obj_ref_mode, 3:] - t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) + as1, bs1, gs1 = np.radians(coord_raw[obj_ref_mode, 3:]) + r_probe = tr.euler_matrix(as1, bs1, gs1, 'rzyx') + t_probe_raw = tr.translation_matrix(coord_raw[obj_ref_mode, :3]) + t_offset_aux = np.linalg.inv(r_s0_raw) @ r_probe @ t_obj_raw + t_offset = np.identity(4) + t_offset[:, -1] = t_offset_aux[:, -1] + t_probe = s0_raw @ t_offset @ np.linalg.inv(s0_raw) @ t_probe_raw + m_probe = tr.concatenate_matrices(t_probe, r_probe) - t_probe_raw[2, -1] = -t_probe_raw[2, -1] + return m_probe - m_img = m_change * t_probe_raw - coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], psi, theta, phi +def object_to_reference(coord_raw, m_probe): + """Compute affine transformation matrix to the reference basis - wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) - - # TODO: Optimize the value of sleep for each tracking device. - sleep(0.175) - - if self._pause_: - return - - -class CoregistrationDynamic(threading.Thread): - """ - Thread to update the coordinates with the fiducial points - co-registration method while the Navigation Button is pressed. - Sleep function in run method is used to avoid blocking GUI and - for better real-time navigation + :param coord_raw: Coordinates returned by the tracking device + :type coord_raw: numpy.ndarray + :param m_probe: Probe coordinates + :type m_probe: numpy.ndarray + :return: 4 x 4 numpy double array + :rtype: numpy.ndarray """ - def __init__(self, coreg_data, nav_id, trck_info): - threading.Thread.__init__(self) - self.coreg_data = coreg_data - self.nav_id = nav_id - self.trck_info = trck_info - self._pause_ = False - self.start() - - def stop(self): - self._pause_ = True - - def run(self): - m_change, obj_ref_mode = self.coreg_data - trck_init, trck_id, trck_mode = self.trck_info - - while self.nav_id: - coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) - - psi, theta, phi = radians(coord_raw[obj_ref_mode, 3:]) - r_probe = tr.euler_matrix(psi, theta, phi, 'rzyx') - t_probe = tr.translation_matrix(coord_raw[obj_ref_mode, :3]) - m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) - - psi_ref, theta_ref, phi_ref = radians(coord_raw[1, 3:]) - r_ref = tr.euler_matrix(psi_ref, theta_ref, phi_ref, 'rzyx') - t_ref = tr.translation_matrix(coord_raw[1, :3]) - m_ref = asmatrix(tr.concatenate_matrices(t_ref, r_ref)) - - m_dyn = m_ref.I * m_probe - m_dyn[2, -1] = -m_dyn[2, -1] - - m_img = m_change * m_dyn - - scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) - - coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ - degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) - - wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) - - # TODO: Optimize the value of sleep for each tracking device. - sleep(0.175) - - if self._pause_: - return - - -class CoregistrationDynamic_old(threading.Thread): - """ - Thread to update the coordinates with the fiducial points - co-registration method while the Navigation Button is pressed. - Sleep function in run method is used to avoid blocking GUI and - for better real-time navigation + a, b, g = np.radians(coord_raw[1, 3:]) + r_ref = tr.euler_matrix(a, b, g, 'rzyx') + t_ref = tr.translation_matrix(coord_raw[1, :3]) + m_ref = tr.concatenate_matrices(t_ref, r_ref) + + m_dyn = np.linalg.inv(m_ref) @ m_probe + return m_dyn + + +def tracker_to_image(m_change, m_probe_ref, r_obj_img, m_obj_raw, s0_dyn): + """Compute affine transformation matrix to the reference basis + + :param m_change: Corregistration transformation obtained from fiducials + :type m_change: numpy.ndarray + :param m_probe_ref: Object or probe in reference coordinate system + :type m_probe_ref: numpy.ndarray + :param r_obj_img: Object coordinate system in image space (3d model) + :type r_obj_img: numpy.ndarray + :param m_obj_raw: Object basis in raw coordinates from tacker + :type m_obj_raw: numpy.ndarray + :param s0_dyn: + :type s0_dyn: numpy.ndarray + :return: 4 x 4 numpy double array + :rtype: numpy.ndarray """ - def __init__(self, bases, nav_id, trck_info): - threading.Thread.__init__(self) - self.bases = bases - self.nav_id = nav_id + m_img = m_change @ m_probe_ref + r_obj = r_obj_img @ np.linalg.inv(m_obj_raw) @ np.linalg.inv(s0_dyn) @ m_probe_ref @ m_obj_raw + m_img[:3, :3] = r_obj[:3, :3] + return m_img + + +def corregistrate_object_dynamic(inp, coord_raw, ref_mode_id): + + m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = inp + + # transform raw marker coordinate to object center + m_probe = object_marker_to_center(coord_raw, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw) + # transform object center to reference marker if specified as dynamic reference + if ref_mode_id: + m_probe_ref = object_to_reference(coord_raw, m_probe) + else: + m_probe_ref = m_probe + # invert y coordinate + m_probe_ref[2, -1] = -m_probe_ref[2, -1] + # corregistrate from tracker to image space + m_img = tracker_to_image(m_change, m_probe_ref, r_obj_img, m_obj_raw, s0_dyn) + # compute rotation angles + _, _, angles, _, _ = tr.decompose_matrix(m_img) + # create output coordiante list + coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ + np.degrees(angles[0]), np.degrees(angles[1]), np.degrees(angles[2]) + + return coord, m_img + + +def compute_marker_transformation(coord_raw, obj_ref_mode): + psi, theta, phi = np.radians(coord_raw[obj_ref_mode, 3:]) + r_probe = tr.euler_matrix(psi, theta, phi, 'rzyx') + t_probe = tr.translation_matrix(coord_raw[obj_ref_mode, :3]) + m_probe = tr.concatenate_matrices(t_probe, r_probe) + return m_probe + + +def corregistrate_dynamic(inp, coord_raw, ref_mode_id): + + m_change, obj_ref_mode = inp + + # transform raw marker coordinate to object center + m_probe = compute_marker_transformation(coord_raw, obj_ref_mode) + # transform object center to reference marker if specified as dynamic reference + if ref_mode_id: + m_ref = compute_marker_transformation(coord_raw, 1) + m_probe_ref = np.linalg.inv(m_ref) @ m_probe + else: + m_probe_ref = m_probe + + # invert y coordinate + m_probe_ref[2, -1] = -m_probe_ref[2, -1] + # corregistrate from tracker to image space + m_img = m_change @ m_probe_ref + # compute rotation angles + _, _, angles, _, _ = tr.decompose_matrix(m_img) + # create output coordiante list + coord = m_img[0, -1], m_img[1, -1], m_img[2, -1],\ + np.degrees(angles[0]), np.degrees(angles[1]), np.degrees(angles[2]) + + return coord, m_img + + +class CoordinateCorregistrate(threading.Thread): + def __init__(self, ref_mode_id, trck_info, coreg_data, coord_queue, view_tracts, coord_tracts_queue, event, sle): + threading.Thread.__init__(self, name='CoordCoregObject') + self.ref_mode_id = ref_mode_id self.trck_info = trck_info - self._pause_ = False - self.start() - - def stop(self): - self._pause_ = True - - def run(self): - m_inv = self.bases[0] - n = self.bases[1] - q1 = self.bases[2] - q2 = self.bases[3] - trck_init = self.trck_info[0] - trck_id = self.trck_info[1] - trck_mode = self.trck_info[2] - - while self.nav_id: - # trck_coord, probe, reference = dco.GetCoordinates(trck_init, trck_id, trck_mode) - coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) - - trck_coord = dco.dynamic_reference(coord_raw[0, :], coord_raw[1, :]) - - trck_xyz = mat([[trck_coord[0]], [trck_coord[1]], [trck_coord[2]]]) - img = q1 + (m_inv * n) * (trck_xyz - q2) - - coord = (float(img[0]), float(img[1]), float(img[2]), trck_coord[3], - trck_coord[4], trck_coord[5]) - angles = coord_raw[0, 3:6] - - # Tried several combinations and different locations to send the messages, - # however only this one does not block the GUI during navigation. - wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=None, position=coord) - wx.CallAfter(Publisher.sendMessage, 'Set camera in volume', coord) - wx.CallAfter(Publisher.sendMessage, 'Update tracker angles', angles) - - # TODO: Optimize the value of sleep for each tracking device. - # Debug tracker is not working with 0.175 so changed to 0.2 - # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. - # sleep(.3) - sleep(0.175) - - if self._pause_: - return - - -class CoregistrationObjectStatic(threading.Thread): - """ - Thread to update the coordinates with the fiducial points - co-registration method while the Navigation Button is pressed. - Sleep function in run method is used to avoid blocking GUI and - for better real-time navigation - """ - - def __init__(self, coreg_data, nav_id, trck_info): - threading.Thread.__init__(self) self.coreg_data = coreg_data - self.nav_id = nav_id - self.trck_info = trck_info - self._pause_ = False - self.start() - - def stop(self): - self._pause_ = True + self.coord_queue = coord_queue + self.view_tracts = view_tracts + self.coord_tracts_queue = coord_tracts_queue + self.event = event + self.sle = sle def run(self): - # m_change = self.coreg_data[0] - # t_obj_raw = self.coreg_data[1] - # s0_raw = self.coreg_data[2] - # r_s0_raw = self.coreg_data[3] - # s0_dyn = self.coreg_data[4] - # m_obj_raw = self.coreg_data[5] - # r_obj_img = self.coreg_data[6] - # obj_ref_mode = self.coreg_data[7] - # - # trck_init = self.trck_info[0] - # trck_id = self.trck_info[1] - # trck_mode = self.trck_info[2] - - m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = self.coreg_data - trck_init, trck_id, trck_mode = self.trck_info - - while self.nav_id: - coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) - - as1, bs1, gs1 = radians(coord_raw[obj_ref_mode, 3:]) - r_probe = asmatrix(tr.euler_matrix(as1, bs1, gs1, 'rzyx')) - t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) - t_offset_aux = r_s0_raw.I * r_probe * t_obj_raw - t_offset = asmatrix(identity(4)) - t_offset[:, -1] = t_offset_aux[:, -1] - t_probe = s0_raw * t_offset * s0_raw.I * t_probe_raw - m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) - - m_probe[2, -1] = -m_probe[2, -1] - - m_img = m_change * m_probe - r_obj = r_obj_img * m_obj_raw.I * s0_dyn.I * m_probe * m_obj_raw - - m_img[:3, :3] = r_obj[:3, :3] - - scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) - - coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ - degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) - - wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) - wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) - - # TODO: Optimize the value of sleep for each tracking device. - sleep(0.175) - - # Debug tracker is not working with 0.175 so changed to 0.2 - # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. - # sleep(.3) - - # partially working for translate and offset, - # but offset is kept always in same axis, have to fix for rotation - # M_dyn = M_reference.I * T_stylus - # M_dyn[2, -1] = -M_dyn[2, -1] - # M_dyn_ch = M_change * M_dyn - # ddd = M_dyn_ch[0, -1], M_dyn_ch[1, -1], M_dyn_ch[2, -1] - # M_dyn_ch[:3, -1] = asmatrix(db.flip_x_m(ddd)).reshape([3, 1]) - # M_final = S0 * M_obj_trans_0 * S0.I * M_dyn_ch - - # this works for static reference object rotation - # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_raw.I * R_stylus * M_obj_rot_raw - # this works for dynamic reference in rotation but not in translation - # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_dyn.I * R_reference.I * R_stylus * M_obj_rot_raw - - if self._pause_: - return - - -class CoregistrationObjectDynamic(threading.Thread): - """ - Thread to update the coordinates with the fiducial points - co-registration method while the Navigation Button is pressed. - Sleep function in run method is used to avoid blocking GUI and - for better real-time navigation - """ - - def __init__(self, coreg_data, nav_id, trck_info): - threading.Thread.__init__(self) - self.coreg_data = coreg_data - self.nav_id = nav_id + trck_info = self.trck_info + coreg_data = self.coreg_data + view_obj = 1 + + trck_init, trck_id, trck_mode = trck_info + # print('CoordCoreg: event {}'.format(self.event.is_set())) + while not self.event.is_set(): + try: + # print(f"Set the coordinate") + coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) + coord, m_img = corregistrate_object_dynamic(coreg_data, coord_raw, self.ref_mode_id) + m_img_flip = m_img.copy() + m_img_flip[1, -1] = -m_img_flip[1, -1] + # self.pipeline.set_message(m_img_flip) + self.coord_queue.put_nowait([coord, m_img, view_obj]) + # print('CoordCoreg: put {}'.format(count)) + # count += 1 + + if self.view_tracts: + self.coord_tracts_queue.put_nowait(m_img_flip) + + # The sleep has to be in both threads + sleep(self.sle) + except queue.Full: + pass + + +class CoordinateCorregistrateNoObject(threading.Thread): + def __init__(self, ref_mode_id, trck_info, coreg_data, coord_queue, view_tracts, coord_tracts_queue, event, sle): + threading.Thread.__init__(self, name='CoordCoregNoObject') + self.ref_mode_id = ref_mode_id self.trck_info = trck_info - self._pause_ = False - self.start() - - def stop(self): - self._pause_ = True + self.coreg_data = coreg_data + self.coord_queue = coord_queue + self.view_tracts = view_tracts + self.coord_tracts_queue = coord_tracts_queue + self.event = event + self.sle = sle def run(self): + trck_info = self.trck_info + coreg_data = self.coreg_data + view_obj = 0 + + trck_init, trck_id, trck_mode = trck_info + # print('CoordCoreg: event {}'.format(self.event.is_set())) + while not self.event.is_set(): + try: + # print(f"Set the coordinate") + coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) + coord, m_img = corregistrate_dynamic(coreg_data, coord_raw, self.ref_mode_id) + # print("Coord: ", coord) + m_img_flip = m_img.copy() + m_img_flip[1, -1] = -m_img_flip[1, -1] + self.coord_queue.put_nowait([coord, m_img, view_obj]) + + if self.view_tracts: + self.coord_tracts_queue.put_nowait(m_img_flip) + + # The sleep has to be in both threads + sleep(self.sle) + except queue.Full: + pass + + +# class CoregistrationStatic(threading.Thread): +# """ +# Thread to update the coordinates with the fiducial points +# co-registration method while the Navigation Button is pressed. +# Sleep function in run method is used to avoid blocking GUI and +# for better real-time navigation +# """ +# +# def __init__(self, coreg_data, nav_id, trck_info): +# threading.Thread.__init__(self) +# self.coreg_data = coreg_data +# self.nav_id = nav_id +# self.trck_info = trck_info +# self._pause_ = False +# self.start() +# +# def stop(self): +# self._pause_ = True +# +# def run(self): +# # m_change = self.coreg_data[0] +# # obj_ref_mode = self.coreg_data[2] +# # +# # trck_init = self.trck_info[0] +# # trck_id = self.trck_info[1] +# # trck_mode = self.trck_info[2] +# +# m_change, obj_ref_mode = self.coreg_data +# trck_init, trck_id, trck_mode = self.trck_info +# +# while self.nav_id: +# coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# +# psi, theta, phi = coord_raw[obj_ref_mode, 3:] +# t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) +# +# t_probe_raw[2, -1] = -t_probe_raw[2, -1] +# +# m_img = m_change * t_probe_raw +# +# coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], psi, theta, phi +# +# wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) +# +# # TODO: Optimize the value of sleep for each tracking device. +# sleep(0.175) +# +# if self._pause_: +# return +# +# +# class CoregistrationDynamic(threading.Thread): +# """ +# Thread to update the coordinates with the fiducial points +# co-registration method while the Navigation Button is pressed. +# Sleep function in run method is used to avoid blocking GUI and +# for better real-time navigation +# """ +# +# def __init__(self, coreg_data, nav_id, trck_info): +# threading.Thread.__init__(self) +# self.coreg_data = coreg_data +# self.nav_id = nav_id +# self.trck_info = trck_info +# # self.tracts_info = tracts_info +# # self.tracts = None +# self._pause_ = False +# # self.start() +# +# def stop(self): +# # self.tracts.stop() +# self._pause_ = True +# +# def run(self): +# m_change, obj_ref_mode = self.coreg_data +# trck_init, trck_id, trck_mode = self.trck_info +# # seed, tracker, affine, affine_vtk = self.tracts_info +# +# while self.nav_id: +# coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# +# psi, theta, phi = radians(coord_raw[obj_ref_mode, 3:]) +# r_probe = tr.euler_matrix(psi, theta, phi, 'rzyx') +# t_probe = tr.translation_matrix(coord_raw[obj_ref_mode, :3]) +# m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) +# +# psi_ref, theta_ref, phi_ref = radians(coord_raw[1, 3:]) +# r_ref = tr.euler_matrix(psi_ref, theta_ref, phi_ref, 'rzyx') +# t_ref = tr.translation_matrix(coord_raw[1, :3]) +# m_ref = asmatrix(tr.concatenate_matrices(t_ref, r_ref)) +# +# m_dyn = m_ref.I * m_probe +# m_dyn[2, -1] = -m_dyn[2, -1] +# +# m_img = m_change * m_dyn +# +# scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) +# +# coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ +# degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) +# +# # pos_world_aux = np.ones([4, 1]) +# # pos_world_aux[:3, -1] = db.flip_x((m_img[0, -1], m_img[1, -1], m_img[2, -1]))[:3] +# # pos_world = np.linalg.inv(affine) @ pos_world_aux +# # seed_aux = pos_world.reshape([1, 4])[0, :3] +# # seed = seed_aux[np.newaxis, :] +# # +# # self.tracts = dtr.compute_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) +# +# # TODO: Optimize the value of sleep for each tracking device. +# sleep(3.175) +# +# if self._pause_: +# return +# +# +# class CoregistrationDynamic_old(threading.Thread): +# """ +# Thread to update the coordinates with the fiducial points +# co-registration method while the Navigation Button is pressed. +# Sleep function in run method is used to avoid blocking GUI and +# for better real-time navigation +# """ +# +# def __init__(self, bases, nav_id, trck_info): +# threading.Thread.__init__(self) +# self.bases = bases +# self.nav_id = nav_id +# self.trck_info = trck_info +# self._pause_ = False +# self.start() +# +# def stop(self): +# self._pause_ = True +# +# def run(self): +# m_inv = self.bases[0] +# n = self.bases[1] +# q1 = self.bases[2] +# q2 = self.bases[3] +# trck_init = self.trck_info[0] +# trck_id = self.trck_info[1] +# trck_mode = self.trck_info[2] +# +# while self.nav_id: +# # trck_coord, probe, reference = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# +# trck_coord = dco.dynamic_reference(coord_raw[0, :], coord_raw[1, :]) +# +# trck_xyz = mat([[trck_coord[0]], [trck_coord[1]], [trck_coord[2]]]) +# img = q1 + (m_inv * n) * (trck_xyz - q2) +# +# coord = (float(img[0]), float(img[1]), float(img[2]), trck_coord[3], +# trck_coord[4], trck_coord[5]) +# angles = coord_raw[0, 3:6] +# +# # Tried several combinations and different locations to send the messages, +# # however only this one does not block the GUI during navigation. +# wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=None, position=coord) +# wx.CallAfter(Publisher.sendMessage, 'Set camera in volume', coord) +# wx.CallAfter(Publisher.sendMessage, 'Update tracker angles', angles) +# +# # TODO: Optimize the value of sleep for each tracking device. +# # Debug tracker is not working with 0.175 so changed to 0.2 +# # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. +# # sleep(.3) +# sleep(0.175) +# +# if self._pause_: +# return +# +# +# class CoregistrationObjectStatic(threading.Thread): +# """ +# Thread to update the coordinates with the fiducial points +# co-registration method while the Navigation Button is pressed. +# Sleep function in run method is used to avoid blocking GUI and +# for better real-time navigation +# """ +# +# def __init__(self, coreg_data, nav_id, trck_info): +# threading.Thread.__init__(self) +# self.coreg_data = coreg_data +# self.nav_id = nav_id +# self.trck_info = trck_info +# self._pause_ = False +# self.start() +# +# def stop(self): +# self._pause_ = True +# +# def run(self): +# # m_change = self.coreg_data[0] +# # t_obj_raw = self.coreg_data[1] +# # s0_raw = self.coreg_data[2] +# # r_s0_raw = self.coreg_data[3] +# # s0_dyn = self.coreg_data[4] +# # m_obj_raw = self.coreg_data[5] +# # r_obj_img = self.coreg_data[6] +# # obj_ref_mode = self.coreg_data[7] +# # +# # trck_init = self.trck_info[0] +# # trck_id = self.trck_info[1] +# # trck_mode = self.trck_info[2] +# +# m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = self.coreg_data +# trck_init, trck_id, trck_mode = self.trck_info +# +# while self.nav_id: +# coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# +# as1, bs1, gs1 = radians(coord_raw[obj_ref_mode, 3:]) +# r_probe = asmatrix(tr.euler_matrix(as1, bs1, gs1, 'rzyx')) +# t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) +# t_offset_aux = r_s0_raw.I * r_probe * t_obj_raw +# t_offset = asmatrix(identity(4)) +# t_offset[:, -1] = t_offset_aux[:, -1] +# t_probe = s0_raw * t_offset * s0_raw.I * t_probe_raw +# m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) +# +# m_probe[2, -1] = -m_probe[2, -1] +# +# m_img = m_change * m_probe +# r_obj = r_obj_img * m_obj_raw.I * s0_dyn.I * m_probe * m_obj_raw +# +# m_img[:3, :3] = r_obj[:3, :3] +# +# scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) +# +# coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ +# degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) +# +# wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) +# wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) +# +# # TODO: Optimize the value of sleep for each tracking device. +# sleep(0.175) +# +# # Debug tracker is not working with 0.175 so changed to 0.2 +# # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. +# # sleep(.3) +# +# # partially working for translate and offset, +# # but offset is kept always in same axis, have to fix for rotation +# # M_dyn = M_reference.I * T_stylus +# # M_dyn[2, -1] = -M_dyn[2, -1] +# # M_dyn_ch = M_change * M_dyn +# # ddd = M_dyn_ch[0, -1], M_dyn_ch[1, -1], M_dyn_ch[2, -1] +# # M_dyn_ch[:3, -1] = asmatrix(db.flip_x_m(ddd)).reshape([3, 1]) +# # M_final = S0 * M_obj_trans_0 * S0.I * M_dyn_ch +# +# # this works for static reference object rotation +# # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_raw.I * R_stylus * M_obj_rot_raw +# # this works for dynamic reference in rotation but not in translation +# # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_dyn.I * R_reference.I * R_stylus * M_obj_rot_raw +# +# if self._pause_: +# return +# +# +# class CoregistrationObjectDynamic(threading.Thread): +# """ +# Thread to update the coordinates with the fiducial points +# co-registration method while the Navigation Button is pressed. +# Sleep function in run method is used to avoid blocking GUI and +# for better real-time navigation +# """ +# +# def __init__(self, coreg_data, nav_id, trck_info, tracts_info): +# threading.Thread.__init__(self) +# self.coreg_data = coreg_data +# self.nav_id = nav_id +# self.trck_info = trck_info +# # self.tracts_info = tracts_info +# # self.tracts = None +# self._pause_ = False +# self.start() +# +# def stop(self): +# # self.tracts.stop() +# self._pause_ = True +# +# def run(self): +# +# m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = self.coreg_data +# trck_init, trck_id, trck_mode = self.trck_info +# # seed, tracker, affine, affine_vtk = self.tracts_info +# +# while self.nav_id: +# coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) +# +# as1, bs1, gs1 = radians(coord_raw[obj_ref_mode, 3:]) +# r_probe = asmatrix(tr.euler_matrix(as1, bs1, gs1, 'rzyx')) +# t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) +# t_offset_aux = r_s0_raw.I * r_probe * t_obj_raw +# t_offset = asmatrix(identity(4)) +# t_offset[:, -1] = t_offset_aux[:, -1] +# t_probe = s0_raw * t_offset * s0_raw.I * t_probe_raw +# m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) +# +# a, b, g = radians(coord_raw[1, 3:]) +# r_ref = tr.euler_matrix(a, b, g, 'rzyx') +# t_ref = tr.translation_matrix(coord_raw[1, :3]) +# m_ref = asmatrix(tr.concatenate_matrices(t_ref, r_ref)) +# +# m_dyn = m_ref.I * m_probe +# m_dyn[2, -1] = -m_dyn[2, -1] +# +# m_img = m_change * m_dyn +# r_obj = r_obj_img * m_obj_raw.I * s0_dyn.I * m_dyn * m_obj_raw +# +# m_img[:3, :3] = r_obj[:3, :3] +# +# scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) +# +# coord = m_img[0, -1], m_img[1, -1], m_img[2, -1],\ +# degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) +# +# # norm_vec = m_img[:3, 2].reshape([1, 3]).tolist() +# # p0 = m_img[:3, -1].reshape([1, 3]).tolist() +# # p2 = [x - 30 * y for x, y in zip(p0[0], norm_vec[0])] +# # m_tract = m_img.copy() +# # m_tract[:3, -1] = np.reshape(np.asarray(p2)[np.newaxis, :], [3, 1]) +# +# # pos_world_aux = np.ones([4, 1]) +# # pos_world_aux[:3, -1] = db.flip_x((p2[0], p2[1], p2[2]))[:3] +# # pos_world = np.linalg.inv(affine) @ pos_world_aux +# # seed_aux = pos_world.reshape([1, 4])[0, :3] +# # seed = seed_aux[np.newaxis, :] +# +# # self.tracts = dtr.compute_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) +# wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) +# +# # TODO: Optimize the value of sleep for each tracking device. +# #sleep(2.175) +# sleep(0.175) +# +# # Debug tracker is not working with 0.175 so changed to 0.2 +# # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. +# # sleep(.3) +# +# # partially working for translate and offset, +# # but offset is kept always in same axis, have to fix for rotation +# # M_dyn = M_reference.I * T_stylus +# # M_dyn[2, -1] = -M_dyn[2, -1] +# # M_dyn_ch = M_change * M_dyn +# # ddd = M_dyn_ch[0, -1], M_dyn_ch[1, -1], M_dyn_ch[2, -1] +# # M_dyn_ch[:3, -1] = asmatrix(db.flip_x_m(ddd)).reshape([3, 1]) +# # M_final = S0 * M_obj_trans_0 * S0.I * M_dyn_ch +# +# # this works for static reference object rotation +# # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_raw.I * R_stylus * M_obj_rot_raw +# # this works for dynamic reference in rotation but not in translation +# # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_dyn.I * R_reference.I * R_stylus * M_obj_rot_raw +# +# if self._pause_: +# return +# +# +# def corregistrate_object(inp, coord_raw): +# m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = inp +# as1, bs1, gs1 = radians(coord_raw[obj_ref_mode, 3:]) +# r_probe = tr.euler_matrix(as1, bs1, gs1, 'rzyx') +# t_probe_raw = tr.translation_matrix(coord_raw[obj_ref_mode, :3]) +# t_offset_aux = np.linalg.inv(r_s0_raw) @ r_probe @ t_obj_raw +# t_offset = identity(4) +# t_offset[:, -1] = t_offset_aux[:, -1] +# t_probe = s0_raw @ t_offset @ np.linalg.inv(s0_raw) @ t_probe_raw +# m_probe = tr.concatenate_matrices(t_probe, r_probe) +# +# a, b, g = radians(coord_raw[1, 3:]) +# r_ref = tr.euler_matrix(a, b, g, 'rzyx') +# t_ref = tr.translation_matrix(coord_raw[1, :3]) +# m_ref = tr.concatenate_matrices(t_ref, r_ref) +# +# m_dyn = np.linalg.inv(m_ref) @ m_probe +# m_dyn[2, -1] = -m_dyn[2, -1] +# +# m_img = m_change @ m_dyn +# r_obj = r_obj_img @ np.linalg.inv(m_obj_raw) @ np.linalg.inv(s0_dyn) @ m_dyn @ m_obj_raw +# +# m_img[:3, :3] = r_obj[:3, :3] +# +# scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) +# +# coord = m_img[0, -1], m_img[1, -1], m_img[2, -1], \ +# degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) +# +# return coord, m_img - m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = self.coreg_data - trck_init, trck_id, trck_mode = self.trck_info - - while self.nav_id: - coord_raw = dco.GetCoordinates(trck_init, trck_id, trck_mode) - - as1, bs1, gs1 = radians(coord_raw[obj_ref_mode, 3:]) - r_probe = asmatrix(tr.euler_matrix(as1, bs1, gs1, 'rzyx')) - t_probe_raw = asmatrix(tr.translation_matrix(coord_raw[obj_ref_mode, :3])) - t_offset_aux = r_s0_raw.I * r_probe * t_obj_raw - t_offset = asmatrix(identity(4)) - t_offset[:, -1] = t_offset_aux[:, -1] - t_probe = s0_raw * t_offset * s0_raw.I * t_probe_raw - m_probe = asmatrix(tr.concatenate_matrices(t_probe, r_probe)) - - a, b, g = radians(coord_raw[1, 3:]) - r_ref = tr.euler_matrix(a, b, g, 'rzyx') - t_ref = tr.translation_matrix(coord_raw[1, :3]) - m_ref = asmatrix(tr.concatenate_matrices(t_ref, r_ref)) - - m_dyn = m_ref.I * m_probe - m_dyn[2, -1] = -m_dyn[2, -1] - - m_img = m_change * m_dyn - r_obj = r_obj_img * m_obj_raw.I * s0_dyn.I * m_dyn * m_obj_raw - - m_img[:3, :3] = r_obj[:3, :3] - - scale, shear, angles, trans, persp = tr.decompose_matrix(m_img) - - coord = m_img[0, -1], m_img[1, -1], m_img[2, -1],\ - degrees(angles[0]), degrees(angles[1]), degrees(angles[2]) - - wx.CallAfter(Publisher.sendMessage, 'Co-registered points', arg=m_img, position=coord) - wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) - - # TODO: Optimize the value of sleep for each tracking device. - sleep(0.175) - - # Debug tracker is not working with 0.175 so changed to 0.2 - # However, 0.2 is too low update frequency ~5 Hz. Need optimization URGENTLY. - # sleep(.3) - - # partially working for translate and offset, - # but offset is kept always in same axis, have to fix for rotation - # M_dyn = M_reference.I * T_stylus - # M_dyn[2, -1] = -M_dyn[2, -1] - # M_dyn_ch = M_change * M_dyn - # ddd = M_dyn_ch[0, -1], M_dyn_ch[1, -1], M_dyn_ch[2, -1] - # M_dyn_ch[:3, -1] = asmatrix(db.flip_x_m(ddd)).reshape([3, 1]) - # M_final = S0 * M_obj_trans_0 * S0.I * M_dyn_ch - - # this works for static reference object rotation - # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_raw.I * R_stylus * M_obj_rot_raw - # this works for dynamic reference in rotation but not in translation - # R_dyn = M_vtk * M_obj_rot_raw.I * S0_rot_dyn.I * R_reference.I * R_stylus * M_obj_rot_raw - - if self._pause_: - return diff --git a/invesalius/data/imagedata_utils.py b/invesalius/data/imagedata_utils.py index e1fbd98..76f8345 100644 --- a/invesalius/data/imagedata_utils.py +++ b/invesalius/data/imagedata_utils.py @@ -37,6 +37,8 @@ from invesalius.data import vtk_utils as vtk_utils import invesalius.reader.bitmap_reader as bitmap_reader import invesalius.utils as utils import invesalius.data.converters as converters +import invesalius.data.slice_ as sl +import invesalius.data.transformations as tr from invesalius import inv_paths @@ -516,3 +518,67 @@ def image_normalize(image, min_=0.0, max_=1.0, output_dtype=np.int16): imin, imax = image.min(), image.max() output[:] = (image - imin) * ((max_ - min_) / (imax - imin)) + min_ return output + + +def world2invspace(affine=None): + """ + Normalize image pixel intensity for int16 gray scale values. + + :param repos: list of translation and rotation [trans_x, trans_y, trans_z, rot_x, rot_y, rot_z] to reposition the + vtk object prior to applying the affine matrix transformation. Note: rotation given in degrees + :param user_matrix: affine matrix from image header, prefered QForm matrix + :return: vtk transform filter for repositioning the polydata and affine matrix to be used as SetUserMatrix in actor + """ + + # remove scaling factor for non-unitary voxel dimensions + scale, shear, angs, trans, persp = tr.decompose_matrix(affine) + affine_noscale = tr.compose_matrix(scale=None, shear=shear, angles=angs, translate=trans, perspective=persp) + # repos_img = [0.] * 6 + # repos_img[1] = -float(shape[1]) + # + # repos_mat = np.identity(4) + # # translation + # repos_mat[:3, -1] = repos_img[:3] + # # rotation (in principle for invesalius space no rotation is needed) + # repos_mat[:3, :3] = tr.euler_matrix(*np.deg2rad(repos_img[3:]), axes='sxyz')[:3, :3] + + # if repos: + # transx, transy, transz, rotx, roty, rotz = repos + # # create a transform that rotates the stl source + # transform = vtk.vtkTransform() + # transform.PostMultiply() + # transform.RotateX(rotx) + # transform.RotateY(roty) + # transform.RotateZ(rotz) + # transform.Translate(transx, transy, transz) + # + # transform_filt = vtk.vtkTransformPolyDataFilter() + # transform_filt.SetTransform(transform) + # transform_filt.Update() + + # assuming vtk default transformation order is PreMultiply, the user matrix is set so: + # 1. Replace the object -> 2. Transform the object to desired position/orientation + # PreMultiplty: M = M*A where M is current transformation and A is applied transformation + # user_matrix = np.linalg.inv(user_matrix) @ repos_mat + + return np.linalg.inv(affine_noscale) + + +def convert_world_to_voxel(xyz, affine): + """ + Convert a coordinate from the world space ((x, y, z); scanner space; millimeters) to the + voxel space ((i, j, k)). This is achieved by multiplying a coordinate by the inverse + of the affine transformation. + More information: https://nipy.org/nibabel/coordinate_systems.html + :param xyz: a list or array of 3 coordinates (x, y, z) in the world coordinates + :param affine: a 4x4 array containing the image affine transformation in homogeneous coordinates + :return: a 1x3 array with the point coordinates in image space (i, j, k) + """ + + # print("xyz: ", xyz, "\naffine", affine) + # convert xyz coordinate to 1x4 homogeneous coordinates array + xyz_homo = np.hstack((xyz, 1.)).reshape([4, 1]) + ijk_homo = np.linalg.inv(affine) @ xyz_homo + ijk = ijk_homo.T[np.newaxis, 0, :3] + + return ijk diff --git a/invesalius/data/record_coords.py b/invesalius/data/record_coords.py index e364cfa..ce29f34 100644 --- a/invesalius/data/record_coords.py +++ b/invesalius/data/record_coords.py @@ -42,7 +42,8 @@ class Record(threading.Thread): self.start() def __bind_events(self): - Publisher.subscribe(self.UpdateCurrentCoords, 'Co-registered points') + # Publisher.subscribe(self.UpdateCurrentCoords, 'Co-registered points') + Publisher.subscribe(self.UpdateCurrentCoords, 'Update cross position') def UpdateCurrentCoords(self, arg, position): self.coord = asarray(position) @@ -50,7 +51,10 @@ class Record(threading.Thread): def stop(self): self._pause_ = True #save coords dialog - filename = dlg.ShowSaveCoordsDialog("coords.csv") + filename = dlg.ShowLoadSaveDialog(message=_(u"Save coords as..."), + wildcard=_("Coordinates files (*.csv)|*.csv"), + style=wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT, + default_filename="coords.csv", save_ext="csv") if filename: savetxt(filename, self.coord_list, delimiter=',', fmt='%.4f', header="time, x, y, z, a, b, g", comments="") diff --git a/invesalius/data/slice_.py b/invesalius/data/slice_.py index a58dee8..ef62bf1 100644 --- a/invesalius/data/slice_.py +++ b/invesalius/data/slice_.py @@ -82,6 +82,9 @@ class Slice(metaclass=utils.Singleton): self.blend_filter = None self.histogram = None self._matrix = None + self._affine = np.identity(4) + self._n_tracts = 0 + self._tracker = None self.aux_matrices = {} self.aux_matrices_colours = {} self.state = const.STATE_DEFAULT @@ -144,6 +147,30 @@ class Slice(metaclass=utils.Singleton): (s * d / 2.0) for (d, s) in zip(self.matrix.shape[::-1], self.spacing) ] + @property + def affine(self): + return self._affine + + @affine.setter + def affine(self, value): + self._affine = value + + @property + def n_tracts(self): + return self._n_tracts + + @n_tracts.setter + def n_tracts(self, value): + self._n_tracts = value + + @property + def tracker(self): + return self._tracker + + @tracker.setter + def tracker(self, value): + self._tracker = value + def __bind_events(self): # General slice control Publisher.subscribe(self.CreateSurfaceFromIndex, "Create surface from index") diff --git a/invesalius/data/styles.py b/invesalius/data/styles.py index 248715b..1b8f7c8 100644 --- a/invesalius/data/styles.py +++ b/invesalius/data/styles.py @@ -47,6 +47,14 @@ from invesalius.data.measures import (CircleDensityMeasure, MeasureData, from invesalius_cy import floodfill +# For tracts +import invesalius.data.tractography as dtr +# import invesalius.project as prj +import invesalius.data.slice_ as sl +import invesalius.data.bases as bases +# from time import sleep +# --- + ORIENTATIONS = { "AXIAL": const.AXIAL, "CORONAL": const.CORONAL, @@ -465,6 +473,18 @@ class CrossInteractorStyle(DefaultInteractorStyle): self.slice_actor = viewer.slice_data.actor self.slice_data = viewer.slice_data + # tracts + # self.seed = [0., 0., 0.] + # slic = sl.Slice() + # self.affine = slic.affine + # self.tracker = slic.tracker + # + # self.affine_vtk = vtk.vtkMatrix4x4() + # for row in range(0, 4): + # for col in range(0, 4): + # self.affine_vtk.SetElement(row, col, self.affine[row, col]) + # --- + self.picker = vtk.vtkWorldPointPicker() self.AddObserver("MouseMoveEvent", self.OnCrossMove) @@ -478,6 +498,7 @@ class CrossInteractorStyle(DefaultInteractorStyle): def CleanUp(self): self.viewer._set_cross_visibility(0) + Publisher.sendMessage('Remove tracts') Publisher.sendMessage('Toggle toolbar item', _id=self.state_code, value=False) @@ -497,12 +518,20 @@ class CrossInteractorStyle(DefaultInteractorStyle): 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) coord = self.viewer.calcultate_scroll_position(px, py) - Publisher.sendMessage('Update cross position', position=(wx, wy, wz)) + + # Tracts + # pos_world_aux = np.ones([4, 1]) + # pos_world_aux[:3, -1] = bases.flip_x((wx, wy, wz))[:3] + # pos_world = np.linalg.inv(self.affine) @ pos_world_aux + # seed_aux = pos_world.reshape([1, 4])[0, :3] + # self.seed = seed_aux[np.newaxis, :] + # print("Check the seed: ", self.seed) + # + + Publisher.sendMessage('Update cross position', arg=None, position=(wx, wy, wz, 0., 0., 0.)) self.ScrollSlice(coord) - Publisher.sendMessage('Set ball reference position', position=(wx, wy, wz)) - Publisher.sendMessage('Co-registered points', arg=None, position=(wx, wy, wz, 0., 0., 0.)) - iren.Render() + # iren.Render() def ScrollSlice(self, coord): if self.orientation == "AXIAL": @@ -522,6 +551,92 @@ class CrossInteractorStyle(DefaultInteractorStyle): index=coord[0]) +class TractsInteractorStyle(CrossInteractorStyle): + """ + Interactor style responsible for tracts visualization. + """ + def __init__(self, viewer): + CrossInteractorStyle.__init__(self, viewer) + + # self.state_code = const.SLICE_STATE_TRACTS + + self.viewer = viewer + # print("Im fucking brilliant!") + self.tracts = None + + # data_dir = b'C:\Users\deoliv1\OneDrive\data\dti' + # FOD_path = b"sub-P0_dwi_FOD.nii" + # full_path = os.path.join(data_dir, FOD_path) + # self.tracker = Trekker.tracker(full_path) + # self.orientation = viewer.orientation + # self.slice_actor = viewer.slice_data.actor + # self.slice_data = viewer.slice_data + + # self.picker = vtk.vtkWorldPointPicker() + + self.AddObserver("MouseMoveEvent", self.OnTractsMove) + self.AddObserver("LeftButtonPressEvent", self.OnTractsMouseClick) + self.AddObserver("LeftButtonReleaseEvent", self.OnTractsReleaseLeftButton) + + # def SetUp(self): + # self.viewer._set_cross_visibility(1) + # Publisher.sendMessage('Toggle toolbar item', + # _id=self.state_code, value=True) + + # def CleanUp(self): + # self.viewer._set_cross_visibility(0) + # Publisher.sendMessage('Toggle toolbar item', + # _id=self.state_code, value=False) + + def OnTractsMove(self, obj, evt): + # The user moved the mouse with left button pressed + if self.left_pressed: + # print("OnTractsMove interactor style") + # iren = obj.GetInteractor() + self.ChangeTracts(True) + + def OnTractsMouseClick(self, obj, evt): + # print("Single mouse click") + # self.tracts = dtr.compute_tracts(self.tracker, self.seed, self.left_pressed) + self.ChangeTracts(True) + + def OnTractsReleaseLeftButton(self, obj, evt): + # time.sleep(3.) + self.tracts.stop() + # self.ChangeCrossPosition(iren) + + def ChangeTracts(self, pressed): + # print("Trying to compute tracts") + self.tracts = dtr.compute_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) + # coord = self.viewer.calcultate_scroll_position(px, py) + # Publisher.sendMessage('Update cross position', position=(wx, wy, wz)) + # # self.ScrollSlice(coord) + # Publisher.sendMessage('Set ball reference position', position=(wx, wy, wz)) + # Publisher.sendMessage('Co-registered points', arg=None, position=(wx, wy, wz, 0., 0., 0.)) + + # iren.Render() + + # def ScrollSlice(self, coord): + # if self.orientation == "AXIAL": + # Publisher.sendMessage(('Set scroll position', 'SAGITAL'), + # index=coord[0]) + # Publisher.sendMessage(('Set scroll position', 'CORONAL'), + # index=coord[1]) + # elif self.orientation == "SAGITAL": + # Publisher.sendMessage(('Set scroll position', 'AXIAL'), + # index=coord[2]) + # Publisher.sendMessage(('Set scroll position', 'CORONAL'), + # index=coord[1]) + # elif self.orientation == "CORONAL": + # Publisher.sendMessage(('Set scroll position', 'AXIAL'), + # index=coord[2]) + # Publisher.sendMessage(('Set scroll position', 'SAGITAL'), + # index=coord[0]) + + class WWWLInteractorStyle(DefaultInteractorStyle): """ Interactor style responsible for Window Level & Width functionality. @@ -2784,6 +2899,7 @@ class Styles: const.SLICE_STATE_SELECT_MASK_PARTS: SelectMaskPartsInteractorStyle, const.SLICE_STATE_FFILL_SEGMENTATION: FloodFillSegmentInteractorStyle, const.SLICE_STATE_CROP_MASK: CropMaskInteractorStyle, + const.SLICE_STATE_TRACTS: TractsInteractorStyle, } @classmethod diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index 80a6ad2..9c454f5 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -383,6 +383,7 @@ class SurfaceManager(): actor = vtk.vtkActor() actor.SetMapper(mapper) + actor.GetProperty().SetBackfaceCulling(1) print("BOunds", actor.GetBounds()) @@ -494,6 +495,7 @@ class SurfaceManager(): # Represent an object (geometry & properties) in the rendered scene actor = vtk.vtkActor() + actor.GetProperty().SetBackfaceCulling(1) actor.SetMapper(mapper) # Set actor colour and transparency @@ -536,6 +538,7 @@ class SurfaceManager(): # Represent an object (geometry & properties) in the rendered scene actor = vtk.vtkActor() + actor.GetProperty().SetBackfaceCulling(1) actor.SetMapper(mapper) del mapper #Create Surface instance diff --git a/invesalius/data/trackers.py b/invesalius/data/trackers.py index b86482e..d2d1f02 100644 --- a/invesalius/data/trackers.py +++ b/invesalius/data/trackers.py @@ -62,6 +62,7 @@ def DefaultTracker(tracker_id): # return tracker initialization variable and type of connection return trck_init, 'wrapper' + def PolarisTracker(tracker_id): from wx import ID_OK trck_init = None @@ -91,6 +92,7 @@ def PolarisTracker(tracker_id): # return tracker initialization variable and type of connection return trck_init, lib_mode + def CameraTracker(tracker_id): trck_init = None try: @@ -105,6 +107,7 @@ def CameraTracker(tracker_id): # return tracker initialization variable and type of connection return trck_init, 'wrapper' + def ClaronTracker(tracker_id): import invesalius.constants as const from invesalius import inv_paths diff --git a/invesalius/data/tractography.py b/invesalius/data/tractography.py new file mode 100644 index 0000000..056d050 --- /dev/null +++ b/invesalius/data/tractography.py @@ -0,0 +1,471 @@ +# -*- coding: utf-8 -*- + +#-------------------------------------------------------------------------- +# Software: InVesalius - Software de Reconstrucao 3D de Imagens Medicas +# Copyright: (C) 2001 Centro de Pesquisas Renato Archer +# Homepage: http://www.softwarepublico.gov.br +# Contact: invesalius@cti.gov.br +# License: GNU - GPL 2 (LICENSE.txt/LICENCA.txt) +#-------------------------------------------------------------------------- +# Este programa e software livre; voce pode redistribui-lo e/ou +# modifica-lo sob os termos da Licenca Publica Geral GNU, conforme +# publicada pela Free Software Foundation; de acordo com a versao 2 +# da Licenca. +# +# Este programa eh distribuido na expectativa de ser util, mas SEM +# QUALQUER GARANTIA; sem mesmo a garantia implicita de +# COMERCIALIZACAO ou de ADEQUACAO A QUALQUER PROPOSITO EM +# PARTICULAR. Consulte a Licenca Publica Geral GNU para obter mais +# detalhes. +#-------------------------------------------------------------------------- + +# Author: Victor Hugo Souza (victorhos-at-hotmail.com) +# Contributions: Dogu Baran Aydogan +# Initial date: 8 May 2020 + +import threading +import time + +import numpy as np +import queue +from pubsub import pub as Publisher +import vtk + +import invesalius.constants as const +import invesalius.data.imagedata_utils as img_utils + +# Nice print for arrays +# np.set_printoptions(precision=2) +# np.set_printoptions(suppress=True) + + +def compute_directions(trk_n): + """Compute direction of a single tract in each point and return as an RGB color + + :param trk_n: nx3 array of doubles (x, y, z) point coordinates composing the tract + :type trk_n: numpy.ndarray + :return: nx3 array of int (x, y, z) RGB colors in the range 0 - 255 + :rtype: numpy.ndarray + """ + + # trk_d = np.diff(trk_n, axis=0, append=2*trk_n[np.newaxis, -1, :]) + trk_d = np.diff(trk_n, axis=0, append=trk_n[np.newaxis, -2, :]) + trk_d[-1, :] *= -1 + # check that linalg norm makes second norm + # https://stackoverflow.com/questions/21030391/how-to-normalize-an-array-in-numpy + direction = 255 * np.absolute((trk_d / np.linalg.norm(trk_d, axis=1)[:, None])) + return direction.astype(int) + + +def compute_tubes(trk, direction): + """Compute and assign colors to a vtkTube for visualization of a single tract + + :param trk: nx3 array of doubles (x, y, z) point coordinates composing the tract + :type trk: numpy.ndarray + :param direction: nx3 array of int (x, y, z) RGB colors in the range 0 - 255 + :type direction: numpy.ndarray + :return: a vtkTubeFilter instance + :rtype: vtkTubeFilter + """ + + numb_points = trk.shape[0] + points = vtk.vtkPoints() + lines = vtk.vtkCellArray() + + colors = vtk.vtkUnsignedCharArray() + colors.SetNumberOfComponents(3) + + k = 0 + lines.InsertNextCell(numb_points) + for j in range(numb_points): + points.InsertNextPoint(trk[j, :]) + colors.InsertNextTuple(direction[j, :]) + lines.InsertCellPoint(k) + k += 1 + + trk_data = vtk.vtkPolyData() + trk_data.SetPoints(points) + trk_data.SetLines(lines) + trk_data.GetPointData().SetScalars(colors) + + # make it a tube + trk_tube = vtk.vtkTubeFilter() + trk_tube.SetRadius(0.5) + trk_tube.SetNumberOfSides(4) + trk_tube.SetInputData(trk_data) + trk_tube.Update() + + return trk_tube + + +def combine_tracts_root(out_list, root, 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 + :rtype: vtkMultiBlockDataSet + """ + + branch = vtk.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): + branch.SetBlock(n, tube.GetOutput()) + + return branch + + +def tracts_computation(trk_list, root, n_tracts): + """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 + :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] + + # 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)] + + root = combine_tracts_root(out_list, root, n_tracts) + + return root + + +def compute_tracts(trekker, position, affine, affine_vtk, n_tracts): + """ Compute tractograms using the Trekker library. + + :param trekker: Trekker library instance + :type trekker: Trekker.T + :param position: 3 double coordinates (x, y, z) in list or array + :type position: list + :param affine: 4 x 4 numpy double array + :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 + """ + + # 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() + # 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', flag=True, root=root, affine_vtk=affine_vtk) + else: + Publisher.sendMessage('Remove tracts') + + +def tracts_computation_branch(trk_list): + """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 + :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) 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 + + +class ComputeTractsThread(threading.Thread): + + def __init__(self, inp, affine_vtk, coord_tracts_queue, tracts_queue, event, sle): + """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. + + 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.coord_tracts_queue = coord_tracts_queue + self.tracts_queue = tracts_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 + # n_threads = n_tracts_total + p_old = np.array([[0., 0., 0.]]) + n_tracts = 0 + + # 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)) + + # 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] + + # 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() + + # print("p_new_shape", coord_offset.shape) + # print("m_img_flip_shape", m_img_flip.shape) + 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: + # 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) + bundle.SetBlock(n_branches, branch) + n_branches += 1 + n_tracts = branch.GetNumberOfBlocks() + + # 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) + if bundle: + bundle.SetBlock(n_branches, branch) + n_tracts += branch.GetNumberOfBlocks() + n_branches += 1 + + else: + bundle = None + + # 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 + # 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, self.affine_vtk, coord_offset)) + # print('ComputeTractsThread: put {}'.format(count)) + + 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) + # 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, 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 + + :param trekker: Trekker instance + :type trekker: Trekker.T + :param params: Dictionary containing the parameters values to set in Trekker. Initial values are in constants.py + :type params: dict + :return: List containing the Trekker instance and number of threads for parallel processing in the computer + :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.writeInterval(params['write_interval']) + trekker.maxLength(params['max_lenth']) + trekker.minLength(params['min_lenth']) + trekker.maxSamplingPerStep(params['max_sampling_step']) + + # 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 = params['numb_threads'] + + trekker.numberOfThreads(n_threads) + # print("Trekker config updated: n_threads, {}; seed_max, {}".format(n_threads, params['seed_max'])) + return trekker, n_threads diff --git a/invesalius/data/trigger.py b/invesalius/data/trigger.py index 081f672..f3272ce 100644 --- a/invesalius/data/trigger.py +++ b/invesalius/data/trigger.py @@ -39,7 +39,7 @@ class Trigger(threading.Thread): try: import serial - self.trigger_init = serial.Serial('COM1', baudrate=9600, timeout=0) + self.trigger_init = serial.Serial('COM5', baudrate=115200, timeout=0) self.COM = True except: @@ -82,3 +82,91 @@ class Trigger(threading.Thread): if self.trigger_init: self.trigger_init.close() return + + +class TriggerNew(threading.Thread): + + def __init__(self, trigger_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='Trigger') + + self.trigger_init = None + self.stylusplh = False + # self.COM = False + self.__bind_events() + try: + import serial + + self.trigger_init = serial.Serial('COM5', baudrate=115200, timeout=0) + # self.COM = True + + except: + #wx.MessageBox(_('Connection with port COM1 failed'), _('Communication error'), wx.OK | wx.ICON_ERROR) + print("Trigger init error: Connection with port COM failed") + # self.COM = False + pass + + # self.coord_queue = coord_queue + self.trigger_queue = trigger_queue + self.event = event + self.sle = sle + + def run(self): + + while not self.event.is_set(): + trigger_on = False + try: + self.trigger_init.write(b'0') + sleep(0.3) + lines = self.trigger_init.readlines() + # Following lines can simulate a trigger in 3 sec repetitions + # sleep(3) + # lines = True + if lines: + trigger_on = True + # wx.CallAfter(Publisher.sendMessage, 'Create marker') + + if self.stylusplh: + trigger_on = True + # wx.CallAfter(Publisher.sendMessage, 'Create marker') + self.stylusplh = False + + self.trigger_queue.put_nowait(trigger_on) + sleep(self.sle) + + except: + print("Trigger not read, error") + pass + + # except queue.Empty: + # pass + # except queue.Full: + # self.coord_queue.task_done() + else: + if self.trigger_init: + self.trigger_init.close() diff --git a/invesalius/data/viewer_slice.py b/invesalius/data/viewer_slice.py index 1ca4ee0..9b06def 100644 --- a/invesalius/data/viewer_slice.py +++ b/invesalius/data/viewer_slice.py @@ -573,7 +573,7 @@ class Viewer(wx.Panel): self.cross.SetFocalPoint((ux, uy, uz)) self.ScrollSlice(coord) - Publisher.sendMessage('Set ball reference position', position=(ux, uy, uz)) + self.interactor.Render() def ScrollSlice(self, coord): if self.orientation == "AXIAL": @@ -817,9 +817,9 @@ class Viewer(wx.Panel): ('Set scroll position', self.orientation)) Publisher.subscribe(self.__update_cross_position, - 'Update cross position') - Publisher.subscribe(self.UpdateSlicesNavigation, - 'Co-registered points') + 'Update cross position') + # Publisher.subscribe(self.UpdateSlicesNavigation, + # 'Co-registered points') ### # Publisher.subscribe(self.ChangeBrushColour, # 'Add mask') @@ -1136,8 +1136,9 @@ class Viewer(wx.Panel): renderer.AddActor(cross_actor) - def __update_cross_position(self, position): - self.cross.SetFocalPoint(position) + def __update_cross_position(self, arg, position): + # self.cross.SetFocalPoint(position[:3]) + self.UpdateSlicesNavigation(None, position) def _set_cross_visibility(self, visibility): self.cross_actor.SetVisibility(visibility) @@ -1296,8 +1297,8 @@ class Viewer(wx.Panel): if self.state == const.SLICE_STATE_CROSS: # Update other slice's cross according to the new focal point from # the actual orientation. - focal_point = self.cross.GetFocalPoint() - Publisher.sendMessage('Update cross position', position=focal_point) + x, y, z = self.cross.GetFocalPoint() + Publisher.sendMessage('Update cross position', arg=None, position=(x, y, z, 0., 0., 0.)) Publisher.sendMessage('Update slice viewer') else: self.interactor.Render() @@ -1509,7 +1510,6 @@ class Viewer(wx.Panel): def UpdateCross(self, coord): self.cross.SetFocalPoint(coord) - Publisher.sendMessage('Set ball reference position', position=self.cross.GetFocalPoint()) Publisher.sendMessage('Co-registered points', arg=None, position=(coord[0], coord[1], coord[2], 0., 0., 0.)) self.OnScrollBar() self.interactor.Render() diff --git a/invesalius/data/viewer_volume.py b/invesalius/data/viewer_volume.py index 27c3bc8..be1c356 100644 --- a/invesalius/data/viewer_volume.py +++ b/invesalius/data/viewer_volume.py @@ -19,9 +19,10 @@ # detalhes. #-------------------------------------------------------------------------- -from math import cos, sin +# from math import cos, sin import os import sys +import time import numpy as np from numpy.core.umath_tests import inner1d @@ -36,6 +37,7 @@ from imageio import imsave import invesalius.constants as const import invesalius.data.bases as bases +import invesalius.data.slice_ as sl import invesalius.data.transformations as tr import invesalius.data.vtk_utils as vtku import invesalius.project as prj @@ -168,6 +170,10 @@ class Viewer(wx.Panel): #self.pTarget = [0., 0., 0.] # self.obj_axes = None + self.x_actor = None + self.y_actor = None + self.z_actor = None + self.mark_actor = None self._mode_cross = False self._to_show_ball = 0 @@ -188,12 +194,31 @@ class Viewer(wx.Panel): self.anglethreshold = const.COIL_ANGLES_THRESHOLD self.distthreshold = const.COIL_COORD_THRESHOLD + # for DTI support tests + # self.ntimes = False + # self._to_show_stream = True + # data_dir = b'C:\Users\deoliv1\OneDrive\data\dti' + # nii_path = b'sub-P0_dwi_FOD.nii' + # trk_path = os.path.join(data_dir, nii_path) + # self.tracker_FOD = Trekker.tracker(trk_path) + # proj = prj.Project() + # self.affine = np.identity(4) + self.actor_tracts = None + self.actor_peel = None + self.seed_offset = const.SEED_OFFSET + # Publisher.sendMessage('Get affine matrix') + + # initialize Trekker parameters + slic = sl.Slice() + affine = slic.affine + self.affine_vtk = vtku.numpy_to_vtkMatrix4x4(affine) + def __bind_events(self): Publisher.subscribe(self.LoadActor, 'Load surface actor into viewer') Publisher.subscribe(self.RemoveActor, 'Remove surface actor from viewer') - Publisher.subscribe(self.OnShowSurface, 'Show surface') + # Publisher.subscribe(self.OnShowSurface, 'Show surface') Publisher.subscribe(self.UpdateRender, 'Render volume viewer') Publisher.subscribe(self.ChangeBackgroundColour, @@ -225,7 +250,8 @@ class Viewer(wx.Panel): Publisher.subscribe(self.LoadSlicePlane, 'Load slice plane') Publisher.subscribe(self.ResetCamClippingRange, 'Reset cam clipping range') - Publisher.subscribe(self.SetVolumeCamera, 'Co-registered points') + # Publisher.subscribe(self.SetVolumeCamera, 'Update cross position') + # Publisher.subscribe(self.SetVolumeCamera, 'Co-registered points') # Publisher.subscribe(self.SetVolumeCamera, 'Set camera in volume') Publisher.subscribe(self.SetVolumeCameraState, 'Update volume camera state') @@ -255,8 +281,8 @@ class Viewer(wx.Panel): Publisher.subscribe(self.RemoveVolume, 'Remove Volume') - Publisher.subscribe(self.SetBallReferencePosition, - 'Set ball reference position') + Publisher.subscribe(self.UpdateCameraBallPosition, + 'Update cross position') Publisher.subscribe(self._check_ball_reference, 'Enable style') Publisher.subscribe(self._uncheck_ball_reference, 'Disable style') @@ -283,11 +309,17 @@ class Viewer(wx.Panel): Publisher.subscribe(self.OnUpdateObjectTargetGuide, 'Update object matrix') Publisher.subscribe(self.OnUpdateTargetCoordinates, 'Update target') Publisher.subscribe(self.OnRemoveTarget, 'Disable or enable coil tracker') - # Publisher.subscribe(self.UpdateObjectTargetView, 'Co-registered points') Publisher.subscribe(self.OnTargetMarkerTransparency, 'Set target transparency') Publisher.subscribe(self.OnUpdateAngleThreshold, 'Update angle threshold') Publisher.subscribe(self.OnUpdateDistThreshold, 'Update dist threshold') + Publisher.subscribe(self.OnUpdateTracts, 'Update tracts') + 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') + def SetStereoMode(self, mode): ren_win = self.interactor.GetRenderWindow() @@ -319,13 +351,23 @@ class Viewer(wx.Panel): def _check_ball_reference(self, style): if style == const.SLICE_STATE_CROSS: self._mode_cross = True - self._check_and_set_ball_visibility() + # self._check_and_set_ball_visibility() + self._ball_ref_visibility = True + # if self._to_show_ball: + if not self.ball_actor: + self.CreateBallReference() + self.interactor.Render() def _uncheck_ball_reference(self, style): if style == const.SLICE_STATE_CROSS: self._mode_cross = False - self.RemoveBallReference() + # self.RemoveBallReference() + self._ball_ref_visibility = False + if self.ball_actor: + self.ren.RemoveActor(self.ball_actor) + self.ball_actor = None + self.interactor.Render() def OnSensors(self, probe_id, ref_id, obj_id=0): @@ -385,12 +427,12 @@ class Viewer(wx.Panel): self.probe = self.ref = self.obj = False self.interactor.Render() - def OnShowSurface(self, index, visibility): - if visibility: - self._to_show_ball += 1 - else: - self._to_show_ball -= 1 - self._check_and_set_ball_visibility() + # def OnShowSurface(self, index, visibility): + # if visibility: + # self._to_show_ball += 1 + # else: + # self._to_show_ball -= 1 + # self._check_and_set_ball_visibility() def OnStartSeed(self): self.seed_points = [] @@ -485,8 +527,8 @@ class Viewer(wx.Panel): if (volumes.GetNumberOfItems()): self.ren.RemoveVolume(volumes.GetLastProp()) self.interactor.Render() - self._to_show_ball -= 1 - self._check_and_set_ball_visibility() + # self._to_show_ball -= 1 + # self._check_and_set_ball_visibility() def RemoveActors(self, actors): "Remove a list of actors" @@ -534,11 +576,11 @@ class Viewer(wx.Panel): Markers created by navigation tools and rendered in volume viewer. """ self.ball_id = ball_id - x, y, z = bases.flip_x(coord) + coord_flip = bases.flip_x_m(coord)[:3, 0] ball_ref = vtk.vtkSphereSource() ball_ref.SetRadius(size) - ball_ref.SetCenter(x, y, z) + ball_ref.SetCenter(coord_flip) mapper = vtk.vtkPolyDataMapper() mapper.SetInputConnection(ball_ref.GetOutputPort()) @@ -557,6 +599,35 @@ class Viewer(wx.Panel): #self.UpdateRender() self.Refresh() + def add_marker(self, coord, color): + """Simplified version for creating a spherical marker in the 3D scene + + :param coord: + :param color: + :return: vtkActor + """ + + # x, y, z = coord + + ball_ref = vtk.vtkSphereSource() + ball_ref.SetRadius(2) + ball_ref.SetCenter(coord) + + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputConnection(ball_ref.GetOutputPort()) + + prop = vtk.vtkProperty() + prop.SetColor(color) + + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.SetProperty(prop) + actor.GetProperty().SetOpacity(.5) + + # ren.AddActor(actor) + + return actor + def HideAllMarkers(self, indexes): ballid = indexes for i in range(0, ballid): @@ -610,7 +681,6 @@ class Viewer(wx.Panel): self.staticballs[index].GetProperty().SetColor(color) self.Refresh() - def OnTargetMarkerTransparency(self, status, index): if status: self.staticballs[index].GetProperty().SetOpacity(1) @@ -1104,7 +1174,6 @@ class Viewer(wx.Panel): cam.SetFocalPoint(cam_focus) cam.SetPosition(cam_pos) - def CreateBallReference(self): """ Red sphere on volume visualization to reference center of @@ -1129,31 +1198,14 @@ class Viewer(wx.Panel): self.ren.AddActor(self.ball_actor) - def ActivateBallReference(self): - self._mode_cross = True - self._ball_ref_visibility = True - if self._to_show_ball: - if not self.ball_actor: - self.CreateBallReference() - - def RemoveBallReference(self): - self._mode_cross = False - self._ball_ref_visibility = False - if self.ball_actor: - self.ren.RemoveActor(self.ball_actor) - self.ball_actor = None + def UpdateCameraBallPosition(self, arg, position): + coord_flip = bases.flip_x_m(position[:3])[:3, 0] + self.ball_actor.SetPosition(coord_flip) + self.SetVolumeCamera(coord_flip) def SetBallReferencePosition(self, position): - if self._to_show_ball: - if not self.ball_actor: - self.ActivateBallReference() - - coord = position - x, y, z = bases.flip_x(coord) - self.ball_actor.SetPosition(x, y, z) - - else: - self.RemoveBallReference() + coord_flip = bases.flip_x_m(position[:3])[:3, 0] + self.ball_actor.SetPosition(coord_flip) def CreateObjectPolyData(self, filename): """ @@ -1225,7 +1277,14 @@ class Viewer(wx.Panel): self.obj_actor.GetProperty().SetOpacity(.4) self.obj_actor.SetVisibility(0) + self.x_actor = self.add_line([0., 0., 0.], [1., 0., 0.], color=[.0, .0, 1.0]) + self.y_actor = self.add_line([0., 0., 0.], [0., 1., 0.], color=[.0, 1.0, .0]) + self.z_actor = self.add_line([0., 0., 0.], [0., 0., 1.], color=[1.0, .0, .0]) + self.ren.AddActor(self.obj_actor) + self.ren.AddActor(self.x_actor) + self.ren.AddActor(self.y_actor) + self.ren.AddActor(self.z_actor) # self.obj_axes = vtk.vtkAxesActor() # self.obj_axes.SetShaftTypeToCylinder() @@ -1236,25 +1295,88 @@ class Viewer(wx.Panel): # self.ren.AddActor(self.obj_axes) - def OnNavigationStatus(self, status): - self.nav_status = status + def add_line(self, p1, p2, color=[0.0, 0.0, 1.0]): + line = vtk.vtkLineSource() + line.SetPoint1(p1) + line.SetPoint2(p2) + + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputConnection(line.GetOutputPort()) + + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetColor(color) + + return actor + + def AddPeeledSurface(self, flag, actor): + if self.actor_peel: + self.ren.RemoveActor(self.actor_peel) + self.actor_peel = None + if flag and actor: + self.ren.AddActor(actor) + self.actor_peel = actor + self.Refresh() + + def OnNavigationStatus(self, nav_status, vis_status): + self.nav_status = nav_status + self.tracts_status = vis_status[1] self.pTarget = self.CenterOfMass() - if self.obj_actor and self.nav_status: - self.obj_actor.SetVisibility(self.obj_state) - if not self.obj_state: - self.Refresh() + + if self.nav_status: + if self.obj_actor: + self.obj_actor.SetVisibility(self.obj_state) + self.x_actor.SetVisibility(self.obj_state) + self.y_actor.SetVisibility(self.obj_state) + self.z_actor.SetVisibility(self.obj_state) + + self.Refresh() + + def UpdateSeedOffset(self, data): + self.seed_offset = data + + def UpdateMarkerOffsetState(self, create=False): + if create: + if not self.mark_actor: + self.mark_actor = self.add_marker([0., 0., 0.], color=[0., 1., 1.]) + self.ren.AddActor(self.mark_actor) + else: + if self.mark_actor: + self.ren.RemoveActor(self.mark_actor) + self.mark_actor = None + self.Refresh() + + def CreateMarkerOffset(self): + self.mark_actor = self.add_marker([0., 0., 0.], color=[0., 1., 1.]) + 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): - m_img[:3, -1] = np.asmatrix(bases.flip_x_m((m_img[0, -1], m_img[1, -1], m_img[2, -1]))).reshape([3, 1]) + # print("Update object orientation") - m_img_vtk = vtk.vtkMatrix4x4() + m_img_flip = m_img.copy() + m_img_flip[1, -1] = -m_img_flip[1, -1] - for row in range(0, 4): - for col in range(0, 4): - m_img_vtk.SetElement(row, col, m_img[row, col]) + # translate coregistered coordinate to display a marker where Trekker seed is computed + # coord_offset = m_img_flip[:3, -1] - self.seed_offset * m_img_flip[:3, 2] + + # print("m_img copy in viewer_vol: {}".format(m_img_copy)) + + # m_img[:3, 0] is from posterior to anterior direction of the coil + # m_img[:3, 1] is from left to right direction of the coil + # m_img[:3, 2] is from bottom to up direction of the coil + + m_img_vtk = vtku.numpy_to_vtkMatrix4x4(m_img_flip) self.obj_actor.SetUserMatrix(m_img_vtk) # self.obj_axes.SetUserMatrix(m_rot_vtk) + self.x_actor.SetUserMatrix(m_img_vtk) + self.y_actor.SetUserMatrix(m_img_vtk) + self.z_actor.SetUserMatrix(m_img_vtk) self.Refresh() @@ -1267,13 +1389,42 @@ class Viewer(wx.Panel): else: if self.obj_actor: self.ren.RemoveActor(self.obj_actor) + self.ren.RemoveActor(self.x_actor) + self.ren.RemoveActor(self.y_actor) + self.ren.RemoveActor(self.z_actor) + self.ren.RemoveActor(self.mark_actor) self.obj_actor = None + self.x_actor = None + self.y_actor = None + self.z_actor = None + self.mark_actor = None self.Refresh() def UpdateShowObjectState(self, state): self.obj_state = state if self.obj_actor and not self.obj_state: self.obj_actor.SetVisibility(self.obj_state) + self.x_actor.SetVisibility(self.obj_state) + self.y_actor.SetVisibility(self.obj_state) + self.z_actor.SetVisibility(self.obj_state) + + self.Refresh() + + def OnUpdateTracts(self, evt=None, flag=None, actor=None, root=None, affine_vtk=None, count=0): + mapper = vtk.vtkCompositePolyDataMapper2() + mapper.SetInputDataObject(root) + + self.actor_tracts = vtk.vtkActor() + self.actor_tracts.SetMapper(mapper) + self.actor_tracts.SetUserMatrix(affine_vtk) + + self.ren.AddActor(self.actor_tracts) + self.Refresh() + + def OnRemoveTracts(self): + if self.actor_tracts: + self.ren.RemoveActor(self.actor_tracts) + self.actor_tracts = None self.Refresh() def __bind_events_wx(self): @@ -1448,10 +1599,12 @@ class Viewer(wx.Panel): def SetVolumeCameraState(self, camera_state): self.camera_state = camera_state - def SetVolumeCamera(self, arg, position): + # def SetVolumeCamera(self, arg, position): + def SetVolumeCamera(self, cam_focus): if self.camera_state: # TODO: exclude dependency on initial focus - cam_focus = np.array(bases.flip_x(position[:3])) + # cam_focus = np.array(bases.flip_x(position[:3])) + # cam_focus = np.array(bases.flip_x(position)) cam = self.ren.GetActiveCamera() if self.initial_focus is None: @@ -1545,16 +1698,16 @@ class Viewer(wx.Panel): def OnShowRaycasting(self): if not self.raycasting_volume: self.raycasting_volume = True - self._to_show_ball += 1 - self._check_and_set_ball_visibility() + # self._to_show_ball += 1 + # self._check_and_set_ball_visibility() if self.on_wl: self.text.Show() def OnHideRaycasting(self): self.raycasting_volume = False self.text.Hide() - self._to_show_ball -= 1 - self._check_and_set_ball_visibility() + # self._to_show_ball -= 1 + # self._check_and_set_ball_visibility() def OnSize(self, evt): self.UpdateRender() @@ -1581,16 +1734,16 @@ class Viewer(wx.Panel): #self.ShowOrientationCube() self.interactor.Render() - self._to_show_ball += 1 - self._check_and_set_ball_visibility() + # self._to_show_ball += 1 + # self._check_and_set_ball_visibility() def RemoveActor(self, actor): utils.debug("RemoveActor") ren = self.ren ren.RemoveActor(actor) self.interactor.Render() - self._to_show_ball -= 1 - self._check_and_set_ball_visibility() + # self._to_show_ball -= 1 + # self._check_and_set_ball_visibility() def RemoveAllActor(self): utils.debug("RemoveAllActor") @@ -1602,7 +1755,8 @@ class Viewer(wx.Panel): def LoadVolume(self, volume, colour, ww, wl): self.raycasting_volume = True - self._to_show_ball += 1 + # self._to_show_ball += 1 + # self._check_and_set_ball_visibility() self.light = self.ren.GetLights().GetNextItem() @@ -1622,15 +1776,14 @@ class Viewer(wx.Panel): self.ren.ResetCamera() self.ren.ResetCameraClippingRange() - self._check_and_set_ball_visibility() self.UpdateRender() def UnloadVolume(self, volume): self.ren.RemoveVolume(volume) del volume self.raycasting_volume = False - self._to_show_ball -= 1 - self._check_and_set_ball_visibility() + # self._to_show_ball -= 1 + # self._check_and_set_ball_visibility() def OnSetViewAngle(self, view): self.SetViewAngle(view) @@ -1777,16 +1930,17 @@ class Viewer(wx.Panel): self.SetViewAngle(const.VOL_ISO) self.repositioned_coronal_plan = 1 - def _check_and_set_ball_visibility(self): - #TODO: When creating Raycasting volume and cross is pressed, it is not - # automatically creating the ball reference. - if self._mode_cross: - if self._to_show_ball > 0 and not self._ball_ref_visibility: - self.ActivateBallReference() - self.interactor.Render() - elif not self._to_show_ball and self._ball_ref_visibility: - self.RemoveBallReference() - self.interactor.Render() + # def _check_and_set_ball_visibility(self): + # #TODO: When creating Raycasting volume and cross is pressed, it is not + # # automatically creating the ball reference. + # # print("mode_cross, show_ball, ball_vis ", self._mode_cross, self._to_show_ball, self._ball_ref_visibility) + # if self._mode_cross: + # if self._to_show_ball > 0 and not self._ball_ref_visibility: + # self.ActivateBallReference() + # self.interactor.Render() + # elif not self._to_show_ball and self._ball_ref_visibility: + # self.RemoveBallReference() + # self.interactor.Render() class SlicePlane: def __init__(self): diff --git a/invesalius/data/vtk_utils.py b/invesalius/data/vtk_utils.py index 87f19d5..ec9c6f0 100644 --- a/invesalius/data/vtk_utils.py +++ b/invesalius/data/vtk_utils.py @@ -287,3 +287,21 @@ class TextZero(object): # font.SetWeight(wx.FONTWEIGHT_BOLD) font.SetSymbolicSize(self.symbolic_syze) canvas.draw_text(self.text, (x, y), font=font) + + +def numpy_to_vtkMatrix4x4(affine): + """ + Convert a numpy 4x4 array to a vtk 4x4 matrix + :param affine: 4x4 array + :return: vtkMatrix4x4 object representing the affine + """ + # test for type and shape of affine matrix + # assert isinstance(affine, np.ndarray) + assert affine.shape == (4, 4) + + affine_vtk = vtk.vtkMatrix4x4() + for row in range(0, 4): + for col in range(0, 4): + affine_vtk.SetElement(row, col, affine[row, col]) + + return affine_vtk diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index 5a8001f..9e4388c 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -387,6 +387,15 @@ def ShowImportOtherFilesDialog(id_type): if id_type == const.ID_NIFTI_IMPORT: dlg.SetMessage(_("Import NIFTi 1 file")) dlg.SetWildcard(WILDCARD_NIFTI) + elif id_type == const.ID_TREKKER_MASK: + dlg.SetMessage(_("Import Trekker mask")) + dlg.SetWildcard(WILDCARD_NIFTI) + elif id_type == const.ID_TREKKER_IMG: + dlg.SetMessage(_("Import Trekker anatomical image")) + dlg.SetWildcard(WILDCARD_NIFTI) + elif id_type == const.ID_TREKKER_FOD: + dlg.SetMessage(_("Import Trekker FOD")) + dlg.SetWildcard(WILDCARD_NIFTI) elif id_type == const.ID_PARREC_IMPORT: dlg.SetMessage(_("Import PAR/REC file")) dlg.SetWildcard(WILDCARD_PARREC) @@ -508,79 +517,11 @@ def ShowSaveAsProjectDialog(default_filename=None): return filename, wildcard == INV_COMPRESSED -# Dialog for neuronavigation markers -def ShowSaveMarkersDialog(default_filename=None): - current_dir = os.path.abspath(".") - dlg = wx.FileDialog(None, - _("Save markers as..."), # title - "", # last used directory - default_filename, - _("Markers files (*.mks)|*.mks"), - wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT) - # dlg.SetFilterIndex(0) # default is VTI - - filename = None - try: - if dlg.ShowModal() == wx.ID_OK: - filename = dlg.GetPath() - ok = 1 - else: - ok = 0 - except(wx._core.PyAssertionError): # TODO: fix win64 - filename = dlg.GetPath() - ok = 1 - - if (ok): - extension = "mks" - if sys.platform != 'win32': - if filename.split(".")[-1] != extension: - filename = filename + "." + extension - - os.chdir(current_dir) - return filename - -def ShowSaveCoordsDialog(default_filename=None): - current_dir = os.path.abspath(".") - dlg = wx.FileDialog(None, - _("Save coords as..."), # title - "", # last used directory - default_filename, - _("Coordinates files (*.csv)|*.csv"), - wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT) - # dlg.SetFilterIndex(0) # default is VTI - - filename = None - try: - if dlg.ShowModal() == wx.ID_OK: - filename = dlg.GetPath() - ok = 1 - else: - ok = 0 - except(wx._core.PyAssertionError): # TODO: fix win64 - filename = dlg.GetPath() - ok = 1 +def ShowLoadSaveDialog(message=_(u"Load File"), current_dir=os.path.abspath("."), style=wx.FD_OPEN | wx.FD_CHANGE_DIR, + wildcard=_("Registration files (*.obr)|*.obr"), default_filename="", save_ext=None): - if (ok): - extension = "csv" - if sys.platform != 'win32': - if filename.split(".")[-1] != extension: - filename = filename + "." + extension - - os.chdir(current_dir) - return filename - - -def ShowLoadMarkersDialog(): - current_dir = os.path.abspath(".") - - dlg = wx.FileDialog(None, message=_("Load markers"), - defaultDir="", - defaultFile="", - wildcard=_("Markers files (*.mks)|*.mks"), - style=wx.FD_OPEN|wx.FD_CHANGE_DIR) - - # inv3 filter is default - dlg.SetFilterIndex(0) + dlg = wx.FileDialog(None, message=message, defaultDir="", defaultFile=default_filename, + wildcard=wildcard, style=style) # Show the dialog and retrieve the user response. If it is the OK response, # process the data. @@ -589,73 +530,25 @@ def ShowLoadMarkersDialog(): if dlg.ShowModal() == wx.ID_OK: # This returns a Python list of files that were selected. filepath = dlg.GetPath() + ok_press = 1 + else: + ok_press = 0 except(wx._core.PyAssertionError): # FIX: win64 filepath = dlg.GetPath() + ok_press = 1 - # Destroy the dialog. Don't do this until you are done with it! - # BAD things can happen otherwise! - dlg.Destroy() - os.chdir(current_dir) - return filepath - - -def ShowSaveRegistrationDialog(default_filename=None): - current_dir = os.path.abspath(".") - dlg = wx.FileDialog(None, - _("Save object registration as..."), # title - "", # last used directory - default_filename, - _("Registration files (*.obr)|*.obr"), - wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT) - # dlg.SetFilterIndex(0) # default is VTI - - filename = None - try: - if dlg.ShowModal() == wx.ID_OK: - filename = dlg.GetPath() - ok = 1 - else: - ok = 0 - except(wx._core.PyAssertionError): # TODO: fix win64 - filename = dlg.GetPath() - ok = 1 - - if (ok): - extension = "obr" + # fix the extension if set different than expected + if save_ext and ok_press: + extension = save_ext if sys.platform != 'win32': - if filename.split(".")[-1] != extension: - filename = filename + "." + extension - - os.chdir(current_dir) - return filename - - -def ShowLoadRegistrationDialog(): - current_dir = os.path.abspath(".") - - dlg = wx.FileDialog(None, message=_("Load object registration"), - defaultDir="", - defaultFile="", - wildcard=_("Registration files (*.obr)|*.obr"), - style=wx.FD_OPEN|wx.FD_CHANGE_DIR) - - # inv3 filter is default - dlg.SetFilterIndex(0) - - # Show the dialog and retrieve the user response. If it is the OK response, - # process the data. - filepath = None - try: - if dlg.ShowModal() == wx.ID_OK: - # This returns a Python list of files that were selected. - filepath = dlg.GetPath() - except(wx._core.PyAssertionError): # FIX: win64 - filepath = dlg.GetPath() + if filepath.split(".")[-1] != extension: + filepath = filepath + "." + extension # Destroy the dialog. Don't do this until you are done with it! # BAD things can happen otherwise! dlg.Destroy() os.chdir(current_dir) + return filepath @@ -920,31 +813,9 @@ def SurfaceSelectionRequiredForDuplication(): # Dialogs for neuronavigation mode -def InvalidFiducials(): - msg = _("Fiducials are invalid. Select all coordinates.") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) - else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() - +# ---------------------------------- -def InvalidObjectRegistration(): - msg = _("Perform coil registration before navigation.") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) - else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() - - -def NavigationTrackerWarning(trck_id, lib_mode): +def ShowNavigationTrackerWarning(trck_id, lib_mode): """ Spatial Tracker connection error """ @@ -976,85 +847,46 @@ def NavigationTrackerWarning(trck_id, lib_mode): dlg.Destroy() -def InvalidMarkersFile(): - msg = _("The TXT file is invalid.") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) - else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() - - -def NoMarkerSelected(): - msg = _("No data selected") +def ShowEnterMarkerID(default): + msg = _("Edit marker ID") if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) + dlg = wx.TextEntryDialog(None, "", msg, defaultValue=default) else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) + dlg = wx.TextEntryDialog(None, msg, "InVesalius 3", value=default) dlg.ShowModal() + result = dlg.GetValue() dlg.Destroy() + return result -def DeleteAllMarkers(): - msg = _("Do you really want to delete all markers?") +def ShowConfirmationDialog(msg=_('Proceed?')): + # msg = _("Do you want to delete all markers?") if sys.platform == 'darwin': dlg = wx.MessageDialog(None, "", msg, wx.OK | wx.CANCEL | wx.ICON_QUESTION) else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", + dlg = wx.MessageDialog(None, msg, "InVesalius 3", wx.OK | wx.CANCEL | wx.ICON_QUESTION) result = dlg.ShowModal() dlg.Destroy() return result -def DeleteTarget(): - msg = _("Target deleted") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) - else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() -def NewTarget(): - msg = _("New target selected") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) - else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() +def ShowColorDialog(color_current): + cdata = wx.ColourData() + cdata.SetColour(wx.Colour(color_current)) + dlg = wx.ColourDialog(None, data=cdata) + dlg.GetColourData().SetChooseFull(True) -def InvalidTargetID(): - msg = _("Sorry, you cannot use 'TARGET' ID") - if sys.platform == 'darwin': - dlg = wx.MessageDialog(None, "", msg, - wx.ICON_INFORMATION | wx.OK) + if dlg.ShowModal() == wx.ID_OK: + color_new = dlg.GetColourData().GetColour().Get(includeAlpha=False) else: - dlg = wx.MessageDialog(None, msg, "InVesalius 3 - Neuronavigator", - wx.ICON_INFORMATION | wx.OK) - dlg.ShowModal() - dlg.Destroy() + color_new = None -def EnterMarkerID(default): - msg = _("Edit marker ID") - if sys.platform == 'darwin': - dlg = wx.TextEntryDialog(None, "", msg, defaultValue=default) - else: - dlg = wx.TextEntryDialog(None, msg, "InVesalius 3", value=default) - dlg.ShowModal() - result = dlg.GetValue() dlg.Destroy() - return result + return color_new + +# ---------------------------------- class NewMask(wx.Dialog): @@ -3569,7 +3401,10 @@ class ObjectCalibrationDialog(wx.Dialog): else: coord = coord_raw[0, :] else: - NavigationTrackerWarning(0, 'choose') + ShowNavigationTrackerWarning(0, 'choose') + + if btn_id == 3: + coord = np.zeros([6,]) # Update text controls with tracker coordinates if coord is not None or np.sum(coord) != 0.0: @@ -3582,7 +3417,7 @@ class ObjectCalibrationDialog(wx.Dialog): self.ball_actors[btn_id].GetProperty().SetColor(0.0, 1.0, 0.0) self.Refresh() else: - NavigationTrackerWarning(0, 'choose') + ShowNavigationTrackerWarning(0, 'choose') def OnChoiceRefMode(self, evt): # When ref mode is changed the tracker coordinates are set to nan diff --git a/invesalius/gui/frame.py b/invesalius/gui/frame.py index 89ef8eb..003e9f0 100644 --- a/invesalius/gui/frame.py +++ b/invesalius/gui/frame.py @@ -1185,12 +1185,13 @@ class MenuBar(wx.MenuBar): else: self.FindItemById(const.ID_GOTO_COORD).Enable(False) - def OnEnableNavigation(self, status): + def OnEnableNavigation(self, nav_status, vis_status): """ Disable mode menu when navigation is on. - :param status: Navigation status + :param nav_status: Navigation status + :param vis_status: Status of items visualization during navigation """ - value = status + value = nav_status if value: self.FindItemById(const.ID_MODE_NAVIGATION).Enable(False) else: diff --git a/invesalius/gui/task_navigator.py b/invesalius/gui/task_navigator.py index 0a7d5a6..42b9449 100644 --- a/invesalius/gui/task_navigator.py +++ b/invesalius/gui/task_navigator.py @@ -18,10 +18,13 @@ #-------------------------------------------------------------------------- from functools import partial +# import os +import queue import sys -import os +import threading import numpy as np +# import Trekker import wx try: @@ -31,22 +34,23 @@ except ImportError: import wx.lib.hyperlink as hl import wx.lib.foldpanelbar as fpb +import wx.lib.colourselect as csel import wx.lib.masked.numctrl from pubsub import pub as Publisher -import wx.lib.colourselect as csel -import wx.lib.platebtn as pbtn - -from math import cos, sin, pi from time import sleep -import invesalius.data.transformations as tr import invesalius.constants as const import invesalius.data.bases as db +# import invesalius.data.brainmesh_handler as brain import invesalius.data.coordinates as dco import invesalius.data.coregistration as dcr +import invesalius.data.slice_ as sl import invesalius.data.trackers as dt +import invesalius.data.tractography as dti +import invesalius.data.transformations as tr import invesalius.data.trigger as trig import invesalius.data.record_coords as rec +import invesalius.data.vtk_utils as vtk_utils import invesalius.gui.dialogs as dlg from invesalius import utils @@ -136,7 +140,7 @@ class InnerFoldPanel(wx.Panel): # Study this. fold_panel = fpb.FoldPanelBar(self, -1, wx.DefaultPosition, - (10, 290), 0, fpb.FPB_SINGLE_FOLD) + (10, 310), 0, fpb.FPB_SINGLE_FOLD) # Fold panel style style = fpb.CaptionBarStyle() style.SetCaptionStyle(fpb.CAPTIONBAR_GRADIENT_V) @@ -161,15 +165,23 @@ class InnerFoldPanel(wx.Panel): leftSpacing=0, rightSpacing=0) # Fold 3 - Markers panel - item = fold_panel.AddFoldPanel(_("Extra tools"), collapsed=True) + item = fold_panel.AddFoldPanel(_("Markers"), collapsed=True) mtw = MarkersPanel(item) fold_panel.ApplyCaptionStyle(item, style) fold_panel.AddFoldPanelWindow(item, mtw, spacing= 0, leftSpacing=0, rightSpacing=0) - # Fold 4 - DBS + # Fold 4 - Tractography panel + item = fold_panel.AddFoldPanel(_("Tractography"), collapsed=True) + otw = TractographyPanel(item) + fold_panel.ApplyCaptionStyle(item, style) + fold_panel.AddFoldPanelWindow(item, otw, spacing=0, + leftSpacing=0, rightSpacing=0) + item.Hide() + + # Fold 5 - DBS self.dbs_item = fold_panel.AddFoldPanel(_("Deep Brain Stimulation"), collapsed=True) dtw = DbsPanel(self.dbs_item) #Atribuir nova var, criar panel @@ -239,8 +251,8 @@ class InnerFoldPanel(wx.Panel): def OnHideDbs(self): self.dbs_item.Hide() - def OnCheckStatus(self, status): - if status: + def OnCheckStatus(self, nav_status, vis_status): + if nav_status: self.checktrigger.Enable(False) self.checkobj.Enable(False) else: @@ -295,6 +307,23 @@ class NeuronavigationPanel(wx.Panel): self.obj_reg = None self.obj_reg_status = False self.track_obj = False + self.event = threading.Event() + + self.coord_queue = QueueCustom(maxsize=1) + # self.visualization_queue = QueueCustom(maxsize=1) + self.trigger_queue = QueueCustom(maxsize=1) + self.coord_tracts_queue = QueueCustom(maxsize=1) + self.tracts_queue = QueueCustom(maxsize=1) + + # Tractography parameters + self.trk_inp = None + self.trekker = None + self.n_threads = None + self.view_tracts = False + self.n_tracts = const.N_TRACTS + self.seed_offset = const.SEED_OFFSET + self.seed_radius = const.SEED_RADIUS + self.sleep_nav = const.SLEEP_NAVIGATION self.tracker_id = const.DEFAULT_TRACKER self.ref_mode_id = const.DEFAULT_REF_MODE @@ -405,10 +434,17 @@ class NeuronavigationPanel(wx.Panel): Publisher.subscribe(self.LoadImageFiducials, 'Load image fiducials') Publisher.subscribe(self.UpdateTriggerState, 'Update trigger state') Publisher.subscribe(self.UpdateTrackObjectState, 'Update track object state') - Publisher.subscribe(self.UpdateImageCoordinates, 'Set ball reference position') + Publisher.subscribe(self.UpdateImageCoordinates, 'Update cross position') Publisher.subscribe(self.OnDisconnectTracker, 'Disconnect tracker') Publisher.subscribe(self.UpdateObjectRegistration, 'Update object registration') Publisher.subscribe(self.OnCloseProject, 'Close project data') + Publisher.subscribe(self.UpdateTrekkerObject, 'Update Trekker object') + Publisher.subscribe(self.UpdateNumTracts, 'Update number of tracts') + Publisher.subscribe(self.UpdateSeedOffset, 'Update seed offset') + Publisher.subscribe(self.UpdateSeedRadius, 'Update seed radius') + Publisher.subscribe(self.UpdateSleep, 'Update sleep') + Publisher.subscribe(self.UpdateNumberThreads, 'Update number of threads') + Publisher.subscribe(self.UpdateTractsVisualization, 'Update tracts visualization') def LoadImageFiducials(self, marker_id, coord): for n in const.BTNS_IMG_MKS: @@ -420,7 +456,37 @@ class NeuronavigationPanel(wx.Panel): for m in [0, 1, 2]: self.numctrls_coord[btn_id][m].SetValue(coord[m]) - def UpdateImageCoordinates(self, position): + def UpdateFRE(self, fre): + # TODO: Exhibit FRE in a warning dialog and only starts navigation after user clicks ok + self.txtctrl_fre.SetValue(str(round(fre, 2))) + if fre <= 3: + self.txtctrl_fre.SetBackgroundColour('GREEN') + else: + self.txtctrl_fre.SetBackgroundColour('RED') + + def UpdateTrekkerObject(self, data): + # self.trk_inp = data + self.trekker = data + + def UpdateNumTracts(self, data): + self.n_tracts = data + + def UpdateSeedOffset(self, data): + self.seed_offset = data + + def UpdateSeedRadius(self, data): + self.seed_radius = data + + def UpdateSleep(self, data): + self.sleep_nav = data + + def UpdateNumberThreads(self, data): + self.n_threads = data + + def UpdateTractsVisualization(self, data): + self.view_tracts = data + + def UpdateImageCoordinates(self, arg, position): # TODO: Change from world coordinates to matrix coordinates. They are better for multi software communication. self.current_coord = position for m in [0, 1, 2]: @@ -473,7 +539,7 @@ class NeuronavigationPanel(wx.Panel): label=_("Tracker disconnected successfully")) self.trk_init = dt.TrackerConnection(self.tracker_id, None, 'connect') if not self.trk_init[0]: - dlg.NavigationTrackerWarning(self.tracker_id, self.trk_init[1]) + dlg.ShowNavigationTrackerWarning(self.tracker_id, self.trk_init[1]) ctrl.SetSelection(0) print("Tracker not connected!") else: @@ -488,7 +554,7 @@ class NeuronavigationPanel(wx.Panel): self.trk_init = dt.TrackerConnection(self.tracker_id, trck, 'disconnect') if not self.trk_init[0]: if evt is not False: - dlg.NavigationTrackerWarning(self.tracker_id, 'disconnect') + dlg.ShowNavigationTrackerWarning(self.tracker_id, 'disconnect') self.tracker_id = 0 ctrl.SetSelection(self.tracker_id) Publisher.sendMessage('Update status text in GUI', @@ -507,7 +573,7 @@ class NeuronavigationPanel(wx.Panel): self.tracker_id = choice self.trk_init = dt.TrackerConnection(self.tracker_id, None, 'connect') if not self.trk_init[0]: - dlg.NavigationTrackerWarning(self.tracker_id, self.trk_init[1]) + dlg.ShowNavigationTrackerWarning(self.tracker_id, self.trk_init[1]) self.tracker_id = 0 ctrl.SetSelection(self.tracker_id) @@ -549,7 +615,18 @@ class NeuronavigationPanel(wx.Panel): coord = None if self.trk_init and self.tracker_id: + # if self.tracker_id == const.DEBUGTRACK: + # if btn_id == 3: + # coord1 = np.array([-120., 0., 0., 0., 0., 0.]) + # elif btn_id == 4: + # coord1 = np.array([120., 0., 0., 0., 0., 0.]) + # elif btn_id == 5: + # coord1 = np.array([0., 120., 0., 0., 0., 0.]) + # coord2 = np.zeros([3, 6]) + # coord_raw = np.vstack([coord1, coord2]) + # else: coord_raw = dco.GetCoordinates(self.trk_init, self.tracker_id, self.ref_mode_id) + if self.ref_mode_id: coord = dco.dynamic_reference_m(coord_raw[0, :], coord_raw[1, :]) else: @@ -557,7 +634,7 @@ class NeuronavigationPanel(wx.Panel): coord[2] = -coord[2] else: - dlg.NavigationTrackerWarning(0, 'choose') + dlg.ShowNavigationTrackerWarning(0, 'choose') # Update number controls with tracker coordinates if coord is not None: @@ -569,112 +646,152 @@ class NeuronavigationPanel(wx.Panel): btn_nav = btn[0] choice_trck = btn[1] choice_ref = btn[2] + errors = False + + # initialize jobs list + jobs_list = [] + vis_components = [self.trigger_state, self.view_tracts] + vis_queues = [self.coord_queue, self.trigger_queue, self.tracts_queue] nav_id = btn_nav.GetValue() - if nav_id: + if not nav_id: + self.event.set() + + # print("coord unfinished: {}, queue {}", self.coord_queue.unfinished_tasks, self.coord_queue.qsize()) + # print("coord_tracts unfinished: {}, queue {}", self.coord_tracts_queue.unfinished_tasks, self.coord_tracts_queue.qsize()) + # print("tracts unfinished: {}, queue {}", self.tracts_queue.unfinished_tasks, self.tracts_queue.qsize()) + self.coord_queue.clear() + # self.visualization_queue.clear() + if self.trigger_state: + self.trigger_queue.clear() + if self.view_tracts: + self.coord_tracts_queue.clear() + self.tracts_queue.clear() + + # print("coord after unfinished: {}, queue {}", self.coord_queue.unfinished_tasks, self.coord_queue.qsize()) + # print("coord_tracts after unfinished: {}, queue {}", self.coord_tracts_queue.unfinished_tasks, self.coord_tracts_queue.qsize()) + # print("tracts after unfinished: {}, queue {}", self.tracts_queue.unfinished_tasks, self.tracts_queue.qsize()) + self.coord_queue.join() + # self.visualization_queue.join() + if self.trigger_state: + self.trigger_queue.join() + if self.view_tracts: + self.coord_tracts_queue.join() + self.tracts_queue.join() + + # print("coord join unfinished: {}, queue {}", self.coord_queue.unfinished_tasks, self.coord_queue.qsize()) + # print("vis join unfinished: {}, queue {}", self.visualization_queue.unfinished_tasks, self.visualization_queue.qsize()) + + tooltip = wx.ToolTip(_("Start neuronavigation")) + btn_nav.SetToolTip(tooltip) + + # Enable all navigation buttons + choice_ref.Enable(True) + choice_trck.Enable(True) + for btn_c in self.btns_coord: + btn_c.Enable(True) + + # if self.trigger_state: + # self.trigger.stop() + + Publisher.sendMessage("Navigation status", nav_status=False, vis_status=vis_components) + + else: + if np.isnan(self.fiducials).any(): - dlg.InvalidFiducials() + wx.MessageBox(_("Invalid fiducials, select all coordinates."), _("InVesalius 3")) btn_nav.SetValue(False) - elif not self.trk_init[0]: - dlg.NavigationTrackerWarning(0, 'choose') + elif not self.trk_init[0] or not self.tracker_id: + dlg.ShowNavigationTrackerWarning(0, 'choose') + errors = True else: + if self.event.is_set(): + self.event.clear() + + # prepare GUI for navigation + Publisher.sendMessage("Navigation status", nav_status=True, vis_status=vis_components) + Publisher.sendMessage("Toggle Cross", id=const.SLICE_STATE_CROSS) + Publisher.sendMessage("Hide current mask") tooltip = wx.ToolTip(_("Stop neuronavigation")) btn_nav.SetToolTip(tooltip) - # Disable all navigation buttons + # disable all navigation buttons choice_ref.Enable(False) choice_trck.Enable(False) for btn_c in self.btns_coord: btn_c.Enable(False) - # fids_head_img = np.zeros([3, 3]) - # for ic in range(0, 3): - # fids_head_img[ic, :] = np.asarray(db.flip_x_m(self.fiducials[ic, :])) - # - # m_head_aux, q_head_aux, m_inv_head_aux = db.base_creation(fids_head_img) - # m_head = np.asmatrix(np.identity(4)) - # m_head[:3, :3] = m_head_aux[:3, :3] - - m, q1, minv = db.base_creation_old(self.fiducials[:3, :]) - n, q2, ninv = db.base_creation_old(self.fiducials[3:, :]) - + # fiducials matrix m_change = tr.affine_matrix_from_points(self.fiducials[3:, :].T, self.fiducials[:3, :].T, shear=False, scale=False) - - # coreg_data = [m_change, m_head] - + # initialize spatial tracker parameters tracker_mode = self.trk_init, self.tracker_id, self.ref_mode_id - # FIXME: FRE is taking long to calculate so it updates on GUI delayed to navigation - I think its fixed - # TODO: Exhibit FRE in a warning dialog and only starts navigation after user clicks ok - fre = db.calculate_fre(self.fiducials, minv, n, q1, q2) - - self.txtctrl_fre.SetValue(str(round(fre, 2))) - if fre <= 3: - self.txtctrl_fre.SetBackgroundColour('GREEN') - else: - self.txtctrl_fre.SetBackgroundColour('RED') - - if self.trigger_state: - self.trigger = trig.Trigger(nav_id) - Publisher.sendMessage("Navigation status", status=True) - Publisher.sendMessage("Toggle Cross", id=const.SLICE_STATE_CROSS) - Publisher.sendMessage("Hide current mask") + # compute fiducial registration error (FRE) + # this is the old way to compute the fre, left here to recheck if new works fine. + # fre = db.calculate_fre(self.fiducials, minv, n, q1, q2) + fre = db.calculate_fre_m(self.fiducials) + self.UpdateFRE(fre) if self.track_obj: - if self.obj_reg_status: + # if object tracking is selected + if not self.obj_reg_status: + # check if object registration was performed + wx.MessageBox(_("Perform coil registration before navigation."), _("InVesalius 3")) + errors = True + else: + # if object registration was correctly performed continue with navigation # obj_reg[0] is object 3x3 fiducial matrix and obj_reg[1] is 3x3 orientation matrix obj_fiducials, obj_orients, obj_ref_mode, obj_name = self.obj_reg - if self.trk_init and self.tracker_id: - - coreg_data = [m_change, obj_ref_mode] - - if self.ref_mode_id: - coord_raw = dco.GetCoordinates(self.trk_init, self.tracker_id, self.ref_mode_id) - obj_data = db.object_registration(obj_fiducials, obj_orients, coord_raw, m_change) - coreg_data.extend(obj_data) - - self.correg = dcr.CoregistrationObjectDynamic(coreg_data, nav_id, tracker_mode) - else: - coord_raw = np.array([None]) - obj_data = db.object_registration(obj_fiducials, obj_orients, coord_raw, m_change) - coreg_data.extend(obj_data) - - self.correg = dcr.CoregistrationObjectStatic(coreg_data, nav_id, tracker_mode) + coreg_data = [m_change, obj_ref_mode] + if self.ref_mode_id: + coord_raw = dco.GetCoordinates(self.trk_init, self.tracker_id, self.ref_mode_id) else: - dlg.NavigationTrackerWarning(0, 'choose') - - else: - dlg.InvalidObjectRegistration() - - else: - coreg_data = [m_change, 0] - if self.ref_mode_id: - # self.correg = dcr.CoregistrationDynamic_old(bases_coreg, nav_id, tracker_mode) - self.correg = dcr.CoregistrationDynamic(coreg_data, nav_id, tracker_mode) - else: - self.correg = dcr.CoregistrationStatic(coreg_data, nav_id, tracker_mode) + coord_raw = np.array([None]) - else: - tooltip = wx.ToolTip(_("Start neuronavigation")) - btn_nav.SetToolTip(tooltip) - - # Enable all navigation buttons - choice_ref.Enable(True) - choice_trck.Enable(True) - for btn_c in self.btns_coord: - btn_c.Enable(True) - - if self.trigger_state: - self.trigger.stop() + obj_data = db.object_registration(obj_fiducials, obj_orients, coord_raw, m_change) + coreg_data.extend(obj_data) - self.correg.stop() + jobs_list.append(dcr.CoordinateCorregistrate(self.ref_mode_id, tracker_mode, coreg_data, self.coord_queue, + self.view_tracts, self.coord_tracts_queue, + self.event, self.sleep_nav)) - Publisher.sendMessage("Navigation status", status=False) + else: + coreg_data = (m_change, 0) + jobs_list.append(dcr.CoordinateCorregistrateNoObject(self.ref_mode_id, tracker_mode, coreg_data, + self.coord_queue, + self.view_tracts, self.coord_tracts_queue, + self.event, self.sleep_nav)) + + if not errors: + #TODO: Test the trigger thread + if self.trigger_state: + # self.trigger = trig.Trigger(nav_id) + jobs_list.append(trig.TriggerNew(self.trigger_queue, self.event, self.sleep_nav)) + + if self.view_tracts: + # initialize Trekker parameters + slic = sl.Slice() + affine = slic.affine + 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 + # print("Appending the tract computation thread!") + jobs_list.append(dti.ComputeTractsThread(self.trk_inp, affine_vtk, + self.coord_tracts_queue, self.tracts_queue, + self.event, self.sleep_nav)) + + jobs_list.append(UpdateNavigationScene(vis_queues, vis_components, + self.event, self.sleep_nav)) + + for jobs in jobs_list: + # jobs.daemon = True + jobs.start() + # del jobs def ResetImageFiducials(self): for m in range(0, 3): @@ -699,6 +816,8 @@ class NeuronavigationPanel(wx.Panel): Publisher.sendMessage('Update object registration') Publisher.sendMessage('Update track object state', flag=False, obj_name=False) Publisher.sendMessage('Delete all markers') + Publisher.sendMessage("Update marker offset state", create=False) + Publisher.sendMessage("Remove tracts") # TODO: Reset camera initial focus Publisher.sendMessage('Reset cam clipping range') @@ -830,8 +949,7 @@ class ObjectRegistrationPanel(wx.Panel): def UpdateTrackerInit(self, nav_prop): self.nav_prop = nav_prop - def UpdateNavigationStatus(self, status): - nav_status = status + def UpdateNavigationStatus(self, nav_status, vis_status): if nav_status: self.checkrecordcoords.Enable(1) self.checktrack.Enable(0) @@ -898,37 +1016,50 @@ class ObjectRegistrationPanel(wx.Panel): pass else: - dlg.NavigationTrackerWarning(0, 'choose') + dlg.ShowNavigationTrackerWarning(0, 'choose') def OnLinkLoad(self, event=None): - filename = dlg.ShowLoadRegistrationDialog() + filename = dlg.ShowLoadSaveDialog(message=_(u"Load object registration"), + wildcard=_("Registration files (*.obr)|*.obr")) + # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\pilot_20200131' + # coil_path = 'magstim_coil_dell_laptop.obr' + # filename = os.path.join(data_dir, coil_path) - if filename: - data = np.loadtxt(filename, delimiter='\t') - self.obj_fiducials = data[:, :3] - self.obj_orients = data[:, 3:] - - text_file = open(filename, "r") - header = text_file.readline().split('\t') - text_file.close() - - self.obj_name = header[1] - self.obj_ref_mode = int(header[-1]) - - self.checktrack.Enable(1) - Publisher.sendMessage('Update object registration', - data=(self.obj_fiducials, self.obj_orients, self.obj_ref_mode, self.obj_name)) - Publisher.sendMessage('Update status text in GUI', label=_("Ready")) - self.checktrack.SetValue(True) - 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")) + 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:] + + text_file = open(filename, "r") + header = text_file.readline().split('\t') + text_file.close() + + self.obj_name = header[1] + self.obj_ref_mode = int(header[-1]) + + self.checktrack.Enable(1) + self.checktrack.SetValue(True) + Publisher.sendMessage('Update object registration', + data=(self.obj_fiducials, self.obj_orients, self.obj_ref_mode, self.obj_name)) + Publisher.sendMessage('Update status text in GUI', + 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")) + except: + wx.MessageBox(_("Object registration file incompatible."), _("InVesalius 3")) + Publisher.sendMessage('Update status text in GUI', label="") def ShowSaveObjectDialog(self, evt): if np.isnan(self.obj_fiducials).any() or np.isnan(self.obj_orients).any(): wx.MessageBox(_("Digitize all object fiducials before saving"), _("Save error")) else: - filename = dlg.ShowSaveRegistrationDialog("object_registration.obr") + filename = dlg.ShowLoadSaveDialog(message=_(u"Save object registration as..."), + wildcard=_("Registration files (*.obr)|*.obr"), + style=wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT, + default_filename="object_registration.obr", save_ext="obr") if filename: hdr = 'Object' + "\t" + utils.decode(self.obj_name, const.FS_ENCODE) + "\t" + 'Reference' + "\t" + str('%d' % self.obj_ref_mode) data = np.hstack([self.obj_fiducials, self.obj_orients]) @@ -969,8 +1100,8 @@ class MarkersPanel(wx.Panel): self.tgt_flag = self.tgt_index = None self.nav_status = False - self.marker_colour = (1.0, 1.0, 0.) - self.marker_size = 3 + self.marker_colour = const.MARKER_COLOUR + self.marker_size = const.MARKER_SIZE # Change marker size spin_size = wx.SpinCtrl(self, -1, "", size=wx.Size(40, 23)) @@ -1045,7 +1176,8 @@ class MarkersPanel(wx.Panel): self.Update() def __bind_events(self): - Publisher.subscribe(self.UpdateCurrentCoord, 'Co-registered points') + # Publisher.subscribe(self.UpdateCurrentCoord, 'Co-registered points') + Publisher.subscribe(self.UpdateCurrentCoord, 'Update cross position') Publisher.subscribe(self.OnDeleteSingleMarker, 'Delete fiducial marker') Publisher.subscribe(self.OnDeleteAllMarkers, 'Delete all markers') Publisher.subscribe(self.OnCreateMarker, 'Create marker') @@ -1055,8 +1187,8 @@ class MarkersPanel(wx.Panel): self.current_coord = position[:] #self.current_angle = pubsub_evt.data[1][3:] - def UpdateNavigationStatus(self, status): - if not status: + def UpdateNavigationStatus(self, nav_status, vis_status): + if not nav_status: sleep(0.5) #self.current_coord[3:] = 0, 0, 0 self.nav_status = False @@ -1073,6 +1205,9 @@ class MarkersPanel(wx.Panel): menu_id.AppendSeparator() target_menu = menu_id.Append(1, _('Set as target')) menu_id.Bind(wx.EVT_MENU, self.OnMenuSetTarget, target_menu) + # TODO: Create the remove target option so the user can disable the target without removing the marker + # target_menu_rem = menu_id.Append(3, _('Remove target')) + # menu_id.Bind(wx.EVT_MENU, self.OnMenuRemoveTarget, target_menu_rem) target_menu.Enable(True) self.PopupMenu(menu_id) @@ -1089,10 +1224,10 @@ class MarkersPanel(wx.Panel): if evt == 'TARGET': id_label = evt else: - id_label = dlg.EnterMarkerID(self.lc.GetItemText(list_index, 4)) + id_label = dlg.ShowEnterMarkerID(self.lc.GetItemText(list_index, 4)) if id_label == 'TARGET': id_label = '' - dlg.InvalidTargetID() + wx.MessageBox(_("Invalid TARGET ID."), _("InVesalius 3")) self.lc.SetItem(list_index, 4, id_label) # Add the new ID to exported list if len(self.list_coord[list_index]) > 8: @@ -1122,34 +1257,29 @@ class MarkersPanel(wx.Panel): Publisher.sendMessage('Disable or enable coil tracker', status=True) self.OnMenuEditMarkerId('TARGET') self.tgt_flag = True - dlg.NewTarget() + wx.MessageBox(_("New target selected."), _("InVesalius 3")) def OnMenuSetColor(self, evt): index = self.lc.GetFocusedItem() - cdata = wx.ColourData() - cdata.SetColour(wx.Colour(self.list_coord[index][6]*255,self.list_coord[index][7]*255,self.list_coord[index][8]*255)) - dlg = wx.ColourDialog(self, data=cdata) - dlg.GetColourData().SetChooseFull(True) - if dlg.ShowModal() == wx.ID_OK: - self.r, self.g, self.b = dlg.GetColourData().GetColour().Get(includeAlpha=False) - r = float(self.r) / 255.0 - g = float(self.g) / 255.0 - b = float(self.b) / 255.0 - dlg.Destroy() - color = [r,g,b] - - Publisher.sendMessage('Set new color', index=index, color=color) - - self.list_coord[index][6] = r - self.list_coord[index][7] = g - self.list_coord[index][8] = b + + color_current = [self.list_coord[index][n] * 255 for n in range(6, 9)] + + color_new = dlg.ShowColorDialog(color_current=color_current) + + if color_new: + assert len(color_new) == 3 + for n, col in enumerate(color_new): + self.list_coord[index][n+6] = col/255.0 + + Publisher.sendMessage('Set new color', index=index, color=color_new) def OnDeleteAllMarkers(self, evt=None): if self.list_coord: if evt is None: result = wx.ID_OK else: - result = dlg.DeleteAllMarkers() + # result = dlg.DeleteAllMarkers() + result = dlg.ShowConfirmationDialog(msg=_("Remove all markers? Cannot be undone.")) if result == wx.ID_OK: self.list_coord = [] @@ -1162,7 +1292,7 @@ class MarkersPanel(wx.Panel): self.tgt_flag = self.tgt_index = None Publisher.sendMessage('Disable or enable coil tracker', status=False) if not hasattr(evt, 'data'): - dlg.DeleteTarget() + wx.MessageBox(_("Target deleted."), _("InVesalius 3")) def OnDeleteSingleMarker(self, evt=None, marker_id=None): # OnDeleteSingleMarker is used for both pubsub and button click events @@ -1183,14 +1313,16 @@ class MarkersPanel(wx.Panel): else: index = None + #TODO: There are bugs when no marker is selected, test and improve if index: if self.tgt_flag and self.tgt_index == index[0]: self.tgt_flag = self.tgt_index = None Publisher.sendMessage('Disable or enable coil tracker', status=False) - dlg.DeleteTarget() + wx.MessageBox(_("No data selected."), _("InVesalius 3")) + self.DeleteMarker(index) else: - dlg.NoMarkerSelected() + wx.MessageBox(_("Target deleted."), _("InVesalius 3")) def DeleteMarker(self, index): for i in reversed(index): @@ -1213,12 +1345,16 @@ class MarkersPanel(wx.Panel): self.CreateMarker(self.current_coord, self.marker_colour, self.marker_size) def OnLoadMarkers(self, evt): - filepath = dlg.ShowLoadMarkersDialog() + filename = dlg.ShowLoadSaveDialog(message=_(u"Load markers"), + wildcard=_("Markers files (*.mks)|*.mks")) + # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\pilot_20200131' + # marker_path = 'markers.mks' + # filename = os.path.join(data_dir, marker_path) - if filepath: + if filename: try: count_line = self.lc.GetItemCount() - content = [s.rstrip() for s in open(filepath)] + content = [s.rstrip() for s in open(filename)] for data in content: target = None line = [s for s in data.split()] @@ -1254,7 +1390,7 @@ class MarkersPanel(wx.Panel): self.CreateMarker(coord, colour, size, line[7]) count_line += 1 except: - dlg.InvalidMarkersFile() + wx.MessageBox(_("Invalid markers file."), _("InVesalius 3")) def OnMarkersVisibility(self, evt, ctrl): @@ -1266,7 +1402,11 @@ class MarkersPanel(wx.Panel): ctrl.SetLabel('Hide') def OnSaveMarkers(self, evt): - filename = dlg.ShowSaveMarkersDialog("markers.mks") + filename = dlg.ShowLoadSaveDialog(message=_(u"Save markers as..."), + wildcard=_("Marker files (*.mks)|*.mks"), + style=wx.FD_SAVE | wx.FD_OVERWRITE_PROMPT, + default_filename="markers.mks", save_ext="mks") + if filename: if self.list_coord: text_file = open(filename, "w") @@ -1339,3 +1479,459 @@ class DbsPanel(wx.Panel): except AttributeError: default_colour = wx.SystemSettings_GetColour(wx.SYS_COLOUR_MENUBAR) + +class TractographyPanel(wx.Panel): + + def __init__(self, parent): + wx.Panel.__init__(self, parent) + try: + default_colour = wx.SystemSettings.GetColour(wx.SYS_COLOUR_MENUBAR) + except AttributeError: + default_colour = wx.SystemSettings_GetColour(wx.SYS_COLOUR_MENUBAR) + self.SetBackgroundColour(default_colour) + + self.affine = None + self.affine_vtk = None + self.trekker = None + self.n_tracts = const.N_TRACTS + self.peel_depth = const.PEEL_DEPTH + self.view_tracts = False + self.seed_offset = const.SEED_OFFSET + self.seed_radius = const.SEED_RADIUS + self.sleep_nav = const.SLEEP_NAVIGATION + self.brain_opacity = const.BRAIN_OPACITY + self.brain_peel = None + self.brain_actor = None + self.n_peels = const.MAX_PEEL_DEPTH + self.p_old = np.array([[0., 0., 0.]]) + self.tracts_run = None + self.trekker_cfg = const.TREKKER_CONFIG + self.nav_status = False + + self.SetAutoLayout(1) + self.__bind_events() + + # Button for creating new coil + tooltip = wx.ToolTip(_("Load brain visualization")) + btn_mask = wx.Button(self, -1, _("Load brain"), size=wx.Size(65, 23)) + btn_mask.SetToolTip(tooltip) + btn_mask.Enable(1) + btn_mask.Bind(wx.EVT_BUTTON, self.OnLinkBrain) + # self.btn_new = btn_new + + # Button for import config coil file + tooltip = wx.ToolTip(_("Load FOD")) + btn_load = wx.Button(self, -1, _("Load FOD"), size=wx.Size(65, 23)) + btn_load.SetToolTip(tooltip) + btn_load.Enable(1) + btn_load.Bind(wx.EVT_BUTTON, self.OnLinkFOD) + # self.btn_load = btn_load + + # Save button for object registration + tooltip = wx.ToolTip(_(u"Load Trekker configuration parameters")) + btn_load_cfg = wx.Button(self, -1, _(u"Configure"), size=wx.Size(65, 23)) + btn_load_cfg.SetToolTip(tooltip) + btn_load_cfg.Enable(1) + btn_load_cfg.Bind(wx.EVT_BUTTON, self.OnLoadParameters) + # self.btn_load_cfg = btn_load_cfg + + # Create a horizontal sizer to represent button save + line_btns = wx.BoxSizer(wx.HORIZONTAL) + line_btns.Add(btn_load, 1, wx.LEFT | wx.TOP | wx.RIGHT, 4) + line_btns.Add(btn_load_cfg, 1, wx.LEFT | wx.TOP | wx.RIGHT, 4) + line_btns.Add(btn_mask, 1, wx.LEFT | wx.TOP | wx.RIGHT, 4) + + # Change peeling depth + text_peel_depth = wx.StaticText(self, -1, _("Peeling depth (mm):")) + spin_peel_depth = wx.SpinCtrl(self, -1, "", size=wx.Size(50, 23)) + spin_peel_depth.Enable(1) + spin_peel_depth.SetRange(0, const.MAX_PEEL_DEPTH) + spin_peel_depth.SetValue(const.PEEL_DEPTH) + spin_peel_depth.Bind(wx.EVT_TEXT, partial(self.OnSelectPeelingDepth, ctrl=spin_peel_depth)) + spin_peel_depth.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectPeelingDepth, ctrl=spin_peel_depth)) + + # Change number of tracts + text_ntracts = wx.StaticText(self, -1, _("Number tracts:")) + spin_ntracts = wx.SpinCtrl(self, -1, "", size=wx.Size(50, 23)) + spin_ntracts.Enable(1) + spin_ntracts.SetRange(1, 2000) + spin_ntracts.SetValue(const.N_TRACTS) + spin_ntracts.Bind(wx.EVT_TEXT, partial(self.OnSelectNumTracts, ctrl=spin_ntracts)) + spin_ntracts.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectNumTracts, ctrl=spin_ntracts)) + + # Change seed offset for computing tracts + text_offset = wx.StaticText(self, -1, _("Seed offset (mm):")) + spin_offset = wx.SpinCtrlDouble(self, -1, "", size=wx.Size(50, 23), inc = 0.1) + spin_offset.Enable(1) + spin_offset.SetRange(0, 100.0) + spin_offset.SetValue(self.seed_offset) + spin_offset.Bind(wx.EVT_TEXT, partial(self.OnSelectOffset, ctrl=spin_offset)) + spin_offset.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectOffset, ctrl=spin_offset)) + # self.spin_offset = spin_offset + + # Change seed radius for computing tracts + text_radius = wx.StaticText(self, -1, _("Seed radius (mm):")) + spin_radius = wx.SpinCtrlDouble(self, -1, "", size=wx.Size(50, 23), inc=0.1) + spin_radius.Enable(1) + spin_radius.SetRange(0, 100.0) + spin_radius.SetValue(self.seed_radius) + spin_radius.Bind(wx.EVT_TEXT, partial(self.OnSelectRadius, ctrl=spin_radius)) + spin_radius.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectRadius, ctrl=spin_radius)) + # self.spin_radius = spin_radius + + # Change sleep pause between navigation loops + text_sleep = wx.StaticText(self, -1, _("Sleep (s):")) + spin_sleep = wx.SpinCtrlDouble(self, -1, "", size=wx.Size(50, 23), inc=0.01) + spin_sleep.Enable(1) + spin_sleep.SetRange(0.01, 10.0) + spin_sleep.SetValue(self.sleep_nav) + spin_sleep.Bind(wx.EVT_TEXT, partial(self.OnSelectSleep, ctrl=spin_sleep)) + spin_sleep.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectSleep, ctrl=spin_sleep)) + + # Change opacity of brain mask visualization + text_opacity = wx.StaticText(self, -1, _("Brain opacity:")) + spin_opacity = wx.SpinCtrlDouble(self, -1, "", size=wx.Size(50, 23), inc=0.1) + spin_opacity.Enable(0) + spin_opacity.SetRange(0, 1.0) + spin_opacity.SetValue(self.brain_opacity) + spin_opacity.Bind(wx.EVT_TEXT, partial(self.OnSelectOpacity, ctrl=spin_opacity)) + spin_opacity.Bind(wx.EVT_SPINCTRL, partial(self.OnSelectOpacity, ctrl=spin_opacity)) + self.spin_opacity = spin_opacity + + # Create a horizontal sizer to threshold configs + border = 1 + line_peel_depth = wx.BoxSizer(wx.HORIZONTAL) + line_peel_depth.AddMany([(text_peel_depth, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_peel_depth, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + line_ntracts = wx.BoxSizer(wx.HORIZONTAL) + line_ntracts.AddMany([(text_ntracts, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_ntracts, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + line_offset = wx.BoxSizer(wx.HORIZONTAL) + line_offset.AddMany([(text_offset, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_offset, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + line_radius = wx.BoxSizer(wx.HORIZONTAL) + line_radius.AddMany([(text_radius, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_radius, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + line_sleep = wx.BoxSizer(wx.HORIZONTAL) + line_sleep.AddMany([(text_sleep, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_sleep, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + line_opacity = wx.BoxSizer(wx.HORIZONTAL) + line_opacity.AddMany([(text_opacity, 1, wx.EXPAND | wx.GROW | wx.TOP | wx.RIGHT | wx.LEFT, border), + (spin_opacity, 0, wx.ALL | wx.EXPAND | wx.GROW, border)]) + + # Check box to enable tract visualization + checktracts = wx.CheckBox(self, -1, _('Enable tracts')) + checktracts.SetValue(False) + checktracts.Enable(0) + checktracts.Bind(wx.EVT_CHECKBOX, partial(self.OnEnableTracts, ctrl=checktracts)) + self.checktracts = checktracts + + # Check box to enable surface peeling + checkpeeling = wx.CheckBox(self, -1, _('Peel surface')) + checkpeeling.SetValue(False) + checkpeeling.Enable(0) + checkpeeling.Bind(wx.EVT_CHECKBOX, partial(self.OnShowPeeling, ctrl=checkpeeling)) + self.checkpeeling = checkpeeling + + border_last = 1 + line_checks = wx.BoxSizer(wx.HORIZONTAL) + line_checks.Add(checktracts, 0, wx.ALIGN_LEFT | wx.RIGHT | wx.LEFT, border_last) + line_checks.Add(checkpeeling, 0, wx.ALIGN_RIGHT | wx.RIGHT | wx.LEFT, border_last) + + # Add line sizers into main sizer + border = 1 + border_last = 10 + main_sizer = wx.BoxSizer(wx.VERTICAL) + main_sizer.Add(line_btns, 0, wx.BOTTOM | wx.ALIGN_CENTER_HORIZONTAL, border_last) + main_sizer.Add(line_peel_depth, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_ntracts, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_offset, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_radius, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_sleep, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_opacity, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP, border) + main_sizer.Add(line_checks, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT | wx.TOP | wx.BOTTOM, border_last) + main_sizer.Fit(self) + + self.SetSizer(main_sizer) + self.Update() + + def __bind_events(self): + Publisher.subscribe(self.OnCloseProject, 'Close project data') + Publisher.subscribe(self.OnUpdateTracts, 'Update cross position') + Publisher.subscribe(self.UpdateNavigationStatus, 'Navigation status') + + def OnSelectPeelingDepth(self, evt, ctrl): + self.peel_depth = ctrl.GetValue() + if self.checkpeeling.GetValue(): + actor = self.brain_peel.get_actor(self.peel_depth) + Publisher.sendMessage('Update peel', flag=True, actor=actor) + + def OnSelectNumTracts(self, evt, ctrl): + self.n_tracts = ctrl.GetValue() + # self.tract.n_tracts = ctrl.GetValue() + Publisher.sendMessage('Update number of tracts', data=self.n_tracts) + + def OnSelectOffset(self, evt, ctrl): + self.seed_offset = ctrl.GetValue() + # self.tract.seed_offset = ctrl.GetValue() + Publisher.sendMessage('Update seed offset', data=self.seed_offset) + + def OnSelectRadius(self, evt, ctrl): + self.seed_radius = ctrl.GetValue() + # self.tract.seed_offset = ctrl.GetValue() + Publisher.sendMessage('Update seed radius', data=self.seed_radius) + + def OnSelectSleep(self, evt, ctrl): + self.sleep_nav = ctrl.GetValue() + # self.tract.seed_offset = ctrl.GetValue() + Publisher.sendMessage('Update sleep', data=self.sleep_nav) + + def OnSelectOpacity(self, evt, ctrl): + self.brain_actor.GetProperty().SetOpacity(ctrl.GetValue()) + Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor) + + def OnShowPeeling(self, evt, ctrl): + # self.view_peeling = ctrl.GetValue() + if ctrl.GetValue(): + actor = self.brain_peel.get_actor(self.peel_depth) + else: + actor = None + Publisher.sendMessage('Update peel', flag=ctrl.GetValue(), actor=actor) + + def OnEnableTracts(self, evt, ctrl): + self.view_tracts = ctrl.GetValue() + Publisher.sendMessage('Update tracts visualization', data=self.view_tracts) + if not self.view_tracts: + Publisher.sendMessage('Remove tracts') + Publisher.sendMessage("Update marker offset state", create=False) + + def UpdateNavigationStatus(self, nav_status, vis_status): + 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_TREKKER_MASK) + img_path = dlg.ShowImportOtherFilesDialog(const.ID_TREKKER_IMG) + # data_dir = os.environ.get('OneDriveConsumer') + '\\data\\dti' + # mask_file = 'sub-P0_dwi_mask.nii' + # mask_path = os.path.join(data_dir, mask_file) + # img_file = 'sub-P0_T1w_biascorrected.nii' + # img_path = os.path.join(data_dir, img_file) + + if not self.affine_vtk: + slic = sl.Slice() + self.affine = slic.affine + 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.brain_actor.GetProperty().SetOpacity(self.brain_opacity) + Publisher.sendMessage('Update peel', flag=True, actor=self.brain_actor) + self.checkpeeling.Enable(1) + self.checkpeeling.SetValue(True) + self.spin_opacity.Enable(1) + Publisher.sendMessage('Update status text in GUI', label=_("Brain model loaded")) + except: + 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_TREKKER_FOD) + # Juuso + # data_dir = os.environ.get('OneDriveConsumer') + '\\data\\dti' + # FOD_path = 'sub-P0_dwi_FOD.nii' + # Baran + # data_dir = os.environ.get('OneDrive') + r'\data\dti_navigation\baran\pilot_20200131' + # 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) + + # 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")) + + Publisher.sendMessage('End busy cursor') + + def OnLoadParameters(self, event=None): + import json + filename = dlg.ShowLoadSaveDialog(message=_(u"Load Trekker configuration"), + wildcard=_("JSON file (*.json)|*.json")) + try: + # Check if filename exists, read the JSON file and check if all parameters match + # with the required list defined in the constants module + # if a parameter is missing, raise an error + if filename: + with open(filename) as json_file: + self.trekker_cfg = json.load(json_file) + assert all(name in self.trekker_cfg for name in const.TREKKER_CONFIG) + if self.trekker: + self.trekker, n_threads = dti.set_trekker_parameters(self.trekker, self.trekker_cfg) + Publisher.sendMessage('Update Trekker object', data=self.trekker) + Publisher.sendMessage('Update number of threads', data=n_threads) + + Publisher.sendMessage('Update status text in GUI', label=_("Trekker config loaded")) + + except (AssertionError, json.decoder.JSONDecodeError): + # Inform user that file is not compatible + self.trekker_cfg = const.TREKKER_CONFIG + wx.MessageBox(_("File incompatible, using default configuration."), _("InVesalius 3")) + Publisher.sendMessage('Update status text in GUI', label="") + + def OnUpdateTracts(self, arg, position): + """ + Minimal working version of tract computation. Updates when cross sends Pubsub message to update. + Position refers to the coordinates in InVesalius 2D space. To represent the same coordinates in the 3D space, + flip_x the coordinates and multiply the z coordinate by -1. This is all done in the flix_x function. + + :param arg: event for pubsub + :param position: list or array with the x, y, and z coordinates in InVesalius space + """ + # Minimal working version of tract computation + # It updates when cross updates + # pass + if self.view_tracts and not self.nav_status: + # print("Running during navigation") + coord_flip = db.flip_x_m(position[:3])[:3, 0] + dti.compute_tracts(self.trekker, coord_flip, self.affine, self.affine_vtk, + self.n_tracts) + + def OnCloseProject(self): + self.checktracts.SetValue(False) + self.checktracts.Enable(0) + self.checkpeeling.SetValue(False) + self.checkpeeling.Enable(0) + + self.spin_opacity.SetValue(const.BRAIN_OPACITY) + self.spin_opacity.Enable(0) + Publisher.sendMessage('Update peel', flag=False, actor=self.brain_actor) + + self.peel_depth = const.PEEL_DEPTH + self.n_tracts = const.N_TRACTS + + Publisher.sendMessage('Remove tracts') + + +class QueueCustom(queue.Queue): + """ + A custom queue subclass that provides a :meth:`clear` method. + https://stackoverflow.com/questions/6517953/clear-all-items-from-the-queue + Modified to a LIFO Queue type (Last-in-first-out). Seems to make sense for the navigation + threads, as the last added coordinate should be the first to be processed. + In the first tests in a short run, seems to increase the coord queue size considerably, + possibly limiting the queue size is good. + """ + + def clear(self): + """ + Clears all items from the queue. + """ + + with self.mutex: + unfinished = self.unfinished_tasks - len(self.queue) + if unfinished <= 0: + if unfinished < 0: + raise ValueError('task_done() called too many times') + self.all_tasks_done.notify_all() + self.unfinished_tasks = unfinished + self.queue.clear() + self.not_full.notify_all() + + +class UpdateNavigationScene(threading.Thread): + + def __init__(self, vis_queues, vis_components, event, sle): + """Class (threading) to update the navigation scene with all graphical elements. + + Sleep function in run method is used to avoid blocking GUI and more fluent, real-time navigation + + :param affine_vtk: Affine matrix in vtkMatrix4x4 instance to update objects position in 3D scene + :type affine_vtk: vtkMatrix4x4 + :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='UpdateScene') + self.trigger_state, self.view_tracts = vis_components + self.coord_queue, self.trigger_queue, self.tracts_queue = vis_queues + self.sle = sle + self.event = event + + def run(self): + # count = 0 + while not self.event.is_set(): + got_coords = False + try: + coord, m_img, view_obj = self.coord_queue.get_nowait() + got_coords = True + + # print('UpdateScene: get {}'.format(count)) + + # use of CallAfter is mandatory otherwise crashes the wx interface + if self.view_tracts: + bundle, affine_vtk, coord_offset = self.tracts_queue.get_nowait() + wx.CallAfter(Publisher.sendMessage, 'Remove tracts') + wx.CallAfter(Publisher.sendMessage, 'Update tracts', flag=True, root=bundle, + affine_vtk=affine_vtk) + wx.CallAfter(Publisher.sendMessage, 'Update marker offset', coord_offset=coord_offset) + self.tracts_queue.task_done() + + if self.trigger_state: + trigger_on = self.trigger_queue.get_nowait() + if trigger_on: + wx.CallAfter(Publisher.sendMessage, 'Create marker') + self.trigger_queue.task_done() + + # TODO: If using the view_tracts substitute the raw coord from the offset coordinate, so the user + # see the red cross in the position of the offset marker + wx.CallAfter(Publisher.sendMessage, 'Update cross position', arg=m_img, position=coord) + + if view_obj: + wx.CallAfter(Publisher.sendMessage, 'Update object matrix', m_img=m_img, coord=coord) + + self.coord_queue.task_done() + # print('UpdateScene: done {}'.format(count)) + # count += 1 + + sleep(self.sle) + except queue.Empty: + if got_coords: + self.coord_queue.task_done() + + +class InputAttributes(object): + # taken from https://stackoverflow.com/questions/2466191/set-attributes-from-dictionary-in-python + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) diff --git a/invesalius/project.py b/invesalius/project.py index 2f79620..62ea790 100644 --- a/invesalius/project.py +++ b/invesalius/project.py @@ -313,8 +313,9 @@ class Project(metaclass=Singleton): self.spacing = project["spacing"] if project.get("affine", ""): self.affine = project["affine"] + # self.affine = project.get("affine") Publisher.sendMessage('Update affine matrix', - affine=self.affine, status=True) + affine=np.asarray(self.affine).reshape(4, 4), status=True) self.compress = project.get("compress", True) -- libgit2 0.21.2