diff --git a/icons/align.png b/icons/align.png
new file mode 100644
index 000000000..212d09d6c
Binary files /dev/null and b/icons/align.png differ
diff --git a/icons/stylus.png b/icons/stylus.png
new file mode 100644
index 000000000..f9e47622a
Binary files /dev/null and b/icons/stylus.png differ
diff --git a/invesalius/data/bases.py b/invesalius/data/bases.py
index a60f27cb4..2d72ea1a4 100644
--- a/invesalius/data/bases.py
+++ b/invesalius/data/bases.py
@@ -126,9 +126,7 @@ def calculate_fre(fiducials_raw, fiducials, ref_mode_id, m_change, m_icp=None):
 
     dist = np.zeros([3, 1])
     for i in range(0, 6, 2):
-        p_m, _ = dcr.corregistrate_dynamic(
-            (m_change, 0), fiducials_raw[i : i + 2], ref_mode_id, icp
-        )
+        p_m, _ = dcr.corregistrate_probe(m_change, None, fiducials_raw[i : i + 2], ref_mode_id, icp)
         dist[int(i / 2)] = np.sqrt(np.sum(np.power((p_m[:3] - fiducials[int(i / 2), :]), 2)))
 
     return float(np.sqrt(np.sum(dist**2) / 3))
@@ -199,19 +197,16 @@ 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
-
-    # XXX: Some duplicate processing is done here: the Euler angles are calculated once by
-    #      the lines below, and then again in dco.coordinates_to_transformation_matrix.
-    #
-    a, b, g = np.radians(coords[3, 3:])
-    r_s0_raw = tr.euler_matrix(a, b, g, axes="rzyx")
-
     s0_raw = dco.coordinates_to_transformation_matrix(
         position=coords[3, :3],
         orientation=coords[3, 3:],
         axes="rzyx",
     )
 
+    # copy rotation submatrix from s0_raw
+    r_s0_raw = np.eye(4)
+    r_s0_raw[:3, :3] = s0_raw[:3, :3]
+
     # compute change of basis for object fiducials in source frame
     base_obj_raw, q_obj_raw = base_creation(fids_raw[:3, :3])
     r_obj_raw = np.identity(4)
diff --git a/invesalius/data/coregistration.py b/invesalius/data/coregistration.py
index f5594f46d..73e08b91f 100644
--- a/invesalius/data/coregistration.py
+++ b/invesalius/data/coregistration.py
@@ -48,7 +48,6 @@ def object_marker_to_center(coord_raw, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw
     :return: 4 x 4 numpy double array
     :rtype: numpy.ndarray
     """
-
     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])
@@ -171,13 +170,14 @@ def image_to_tracker(m_change, coord_raw, target, icp, obj_data):
     return m_target_in_tracker
 
 
-def corregistrate_object_dynamic(inp, coord_raw, ref_mode_id, icp):
-    m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = inp
+def corregistrate_probe(m_change, r_stylus, coord_raw, ref_mode_id, icp=[None, None]):
+    if r_stylus is None:
+        # utils.debug("STYLUS ORIENTATION NOT DEFINED!")
+        r_stylus = np.eye(3)
 
-    # 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)
+    m_probe = compute_marker_transformation(coord_raw, 0)
 
-    # transform object center to reference marker if specified as dynamic reference
+    # transform probe to reference system if dynamic ref_mode
     if ref_mode_id:
         m_probe_ref = object_to_reference(coord_raw, m_probe)
     else:
@@ -186,71 +186,109 @@ def corregistrate_object_dynamic(inp, coord_raw, ref_mode_id, icp):
     # 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)
+    # translate m_probe_ref from tracker to image space
+    m_img = m_change @ m_probe_ref
     m_img = apply_icp(m_img, icp)
 
+    # Rotate from trk system where stylus points in x-axis to vtk-system where stylus points in y-axis
+    R = tr.euler_matrix(*np.radians([0, 0, -90]), axes="rxyz")[:3, :3]
+
+    # rotate m_probe_ref from tracker to image space
+    r_img = r_stylus @ R @ m_probe_ref[:3, :3] @ np.linalg.inv(R)
+    m_img[:3, :3] = r_img[:3, :3]
+
     # compute rotation angles
-    angles = tr.euler_from_matrix(m_img, axes="sxyz")
+    angles = np.degrees(tr.euler_from_matrix(m_img, axes="sxyz"))
 
     # create output coordinate 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]),
+        angles[0],
+        angles[1],
+        angles[2],
     )
 
     return coord, m_img
 
 
-def compute_marker_transformation(coord_raw, obj_ref_mode):
-    m_probe = dco.coordinates_to_transformation_matrix(
-        position=coord_raw[obj_ref_mode, :3],
-        orientation=coord_raw[obj_ref_mode, 3:],
-        axes="rzyx",
+def corregistrate_object_dynamic(inp, coord_raw, i_obj, icp):
+    """
+    Corregistrate the object at coord_raw[i_obj] in dynamic ref_mode
+    """
+    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, i_obj, t_obj_raw, s0_raw, r_s0_raw)
+
+    # transform object center to reference marker
+    m_probe_ref = object_to_reference(coord_raw, 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)
+    m_img = apply_icp(m_img, icp)
+
+    # compute rotation angles
+    angles = np.degrees(tr.euler_from_matrix(m_img, axes="sxyz"))
+
+    # create output coordinate list
+    coord = (
+        m_img[0, -1],
+        m_img[1, -1],
+        m_img[2, -1],
+        angles[0],
+        angles[1],
+        angles[2],
     )
-    return m_probe
 
+    return coord, m_img
 
-def corregistrate_dynamic(inp, coord_raw, ref_mode_id, icp):
-    m_change, obj_ref_mode = inp
 
-    # transform raw marker coordinate to object center
-    m_probe = compute_marker_transformation(coord_raw, obj_ref_mode)
+def corregistrate_object_static(inp, coord_raw, i_obj, icp):
+    """
+    Corregistrate the object at coord_raw[i_obj] in static ref_mode
+    """
+    m_change, obj_ref_mode, t_obj_raw, s0_raw, r_s0_raw, s0_dyn, m_obj_raw, r_obj_img = inp
 
-    # 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
+    # transform raw marker coordinate to object center
+    m_probe = object_marker_to_center(coord_raw, i_obj, t_obj_raw, s0_raw, r_s0_raw)
 
     # invert y coordinate
-    m_probe_ref[2, -1] = -m_probe_ref[2, -1]
+    m_probe[2, -1] = -m_probe[2, -1]
 
     # corregistrate from tracker to image space
-    m_img = m_change @ m_probe_ref
+    m_img = tracker_to_image(m_change, m_probe, r_obj_img, m_obj_raw, s0_dyn)
     m_img = apply_icp(m_img, icp)
 
     # compute rotation angles
-    angles = tr.euler_from_matrix(m_img, axes="sxyz")
+    angles = np.degrees(tr.euler_from_matrix(m_img, axes="sxyz"))
 
     # create output coordinate 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]),
+        angles[0],
+        angles[1],
+        angles[2],
     )
 
     return coord, m_img
 
 
+def compute_marker_transformation(coord_raw, obj_ref_mode):
+    m_probe = dco.coordinates_to_transformation_matrix(
+        position=coord_raw[obj_ref_mode, :3],
+        orientation=coord_raw[obj_ref_mode, 3:],
+        axes="rzyx",
+    )
+    return m_probe
+
+
 def apply_icp(m_img, icp):
     use_icp, m_icp = icp
     if use_icp:
@@ -295,6 +333,7 @@ def __init__(
         self,
         ref_mode_id,
         tracker,
+        n_coils,
         coreg_data,
         view_tracts,
         queues,
@@ -308,6 +347,7 @@ def __init__(
         threading.Thread.__init__(self, name="CoordCoregObject")
         self.ref_mode_id = ref_mode_id
         self.tracker = tracker
+        self.n_coils = n_coils
         self.coreg_data = coreg_data
         self.coord_queue = queues[0]
         self.view_tracts = view_tracts
@@ -336,9 +376,20 @@ def __init__(
 
     def run(self):
         coreg_data = self.coreg_data
-        view_obj = 1
+        corregistrate_object = (
+            corregistrate_object_dynamic if self.ref_mode_id else corregistrate_object_static
+        )
+
+        # compute n_coils_effective, the no. of coils to actually process:
+        # check how many coords we get from tracker (-2 for probe & head)
+        n_coils_trk = self.tracker.TrackerCoordinates.GetCoordinates()[0].shape[0] - 2
+
+        obj_ref_mode = coreg_data[2]
+        # if obj_ref_mode=0: only coregister one coil
+        # else: process the other (n_coils - 1) coils too
+        # min(obj_ref_mode, 1) = 0 if obj_ref_mode==0 else 1
+        n_coils_effective = min(n_coils_trk, 1 + min(obj_ref_mode, 1) * (self.n_coils - 1))
 
-        # print('CoordCoreg: event {}'.format(self.event.is_set()))
         while not self.event.is_set():
             try:
                 if not self.icp_queue.empty():
@@ -347,12 +398,27 @@ def run(self):
                 if not self.object_at_target_queue.empty():
                     self.target_flag = self.object_at_target_queue.get_nowait()
 
-                # print(f"Set the coordinate")
                 coord_raw, marker_visibilities = self.tracker.TrackerCoordinates.GetCoordinates()
-                coord, m_img = corregistrate_object_dynamic(
-                    coreg_data, coord_raw, self.ref_mode_id, [self.use_icp, self.m_icp]
+
+                # m_change = coreg_data[1], r_stylus = coreg_data[0]
+                coord_probe, m_img_probe = corregistrate_probe(
+                    coreg_data[1], coreg_data[0], coord_raw, self.ref_mode_id
+                )
+                coord_coil, m_img_coil = corregistrate_object(
+                    coreg_data[1:], coord_raw, obj_ref_mode, [self.use_icp, self.m_icp]
                 )
 
+                coords = [coord_probe, coord_coil]
+                m_imgs = [m_img_probe, m_img_coil]
+
+                # the possible other coils are i_obj=3 onwards at coord_raw[i_obj]
+                for i_obj in range(3, n_coils_effective + 2):
+                    coord_coil, m_img_coil = corregistrate_object(
+                        coreg_data[1:], coord_raw, i_obj, [self.use_icp, self.m_icp]
+                    )
+                    coords.append(coord_coil)
+                    m_imgs.append(m_img_coil)
+
                 # XXX: This is not the best place to do the logic related to approaching the target when the
                 #      debug tracker is in use. However, the trackers (including the debug trackers) operate in
                 #      the tracker space where it is hard to make the tracker approach the target in the image space.
@@ -363,77 +429,26 @@ def run(self):
                 #
                 if self.tracker_id == const.DEBUGTRACKAPPROACH and self.target is not None:
                     if self.last_coord is None:
-                        self.last_coord = np.array(coord)
+                        self.last_coord = np.array(coord_coil)
                     else:
                         coord = self.last_coord + (self.target - self.last_coord) * 0.05
+                        coords[1] = coord
                         self.last_coord = coord
 
                     angles = [np.radians(coord[3]), np.radians(coord[4]), np.radians(coord[5])]
                     translate = coord[0:3]
-                    m_img = tr.compose_matrix(angles=angles, translate=translate)
+                    m_imgs[1] = tr.compose_matrix(angles=angles, translate=translate)
 
-                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([coords, marker_visibilities, m_imgs])
 
-                self.coord_queue.put_nowait([coord, marker_visibilities, m_img, view_obj])
-                # print('CoordCoreg: put {}'.format(count))
-                # count += 1
+                coord = coords[1]  # main coil
+                m_img = m_imgs[1]
+                # LUKATODO: should coord = coords[track_coil]
+                # should the stuff below ever be done for stylus, but not coil?
 
-                if self.view_tracts:
-                    self.coord_tracts_queue.put_nowait(m_img_flip)
-                if self.e_field_loaded:
-                    self.efield_queue.put_nowait([m_img, coord])
-                if not self.icp_queue.empty():
-                    self.icp_queue.task_done()
-            except queue.Full:
-                pass
-            # The sleep has to be in both threads
-            sleep(self.sle)
-
-
-class CoordinateCorregistrateNoObject(threading.Thread):
-    def __init__(
-        self, ref_mode_id, tracker, coreg_data, view_tracts, queues, event, sle, icp, e_field_loaded
-    ):
-        threading.Thread.__init__(self, name="CoordCoregNoObject")
-        self.ref_mode_id = ref_mode_id
-        self.tracker = tracker
-        self.coreg_data = coreg_data
-        self.coord_queue = queues[0]
-        self.view_tracts = view_tracts
-        self.coord_tracts_queue = queues[1]
-        self.event = event
-        self.sle = sle
-        self.icp_queue = queues[2]
-        self.use_icp = icp.use_icp
-        self.m_icp = icp.m_icp
-        self.efield_queue = queues[3]
-        self.e_field_loaded = e_field_loaded
-
-    def run(self):
-        coreg_data = self.coreg_data
-        view_obj = 0
-
-        # print('CoordCoreg: event {}'.format(self.event.is_set()))
-        while not self.event.is_set():
-            try:
-                if self.icp_queue.empty():
-                    None
-                else:
-                    self.use_icp, self.m_icp = self.icp_queue.get_nowait()
-                # print(f"Set the coordinate")
-                # print(self.icp, self.m_icp)
-                coord_raw, marker_visibilities = self.tracker.TrackerCoordinates.GetCoordinates()
-                coord, m_img = corregistrate_dynamic(
-                    coreg_data, coord_raw, self.ref_mode_id, [self.use_icp, self.m_icp]
-                )
-                # 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, marker_visibilities, m_img, view_obj])
-
                 if self.view_tracts:
                     self.coord_tracts_queue.put_nowait(m_img_flip)
                 if self.e_field_loaded:
diff --git a/invesalius/data/viewer_volume.py b/invesalius/data/viewer_volume.py
index 4643683cf..9cf5010a5 100644
--- a/invesalius/data/viewer_volume.py
+++ b/invesalius/data/viewer_volume.py
@@ -1,6 +1,3 @@
-# -*- coding: utf-8 -*-
-import math
-
 # --------------------------------------------------------------------------
 # Software:     InVesalius - Software de Reconstrucao 3D de Imagens Medicas
 # Copyright:    (C) 2001  Centro de Pesquisas Renato Archer
@@ -22,7 +19,6 @@
 # from math import cos, sin
 import os
 import queue
-import random
 import sys
 
 import numpy as np
@@ -33,31 +29,24 @@
 from vtkmodules.vtkCommonColor import vtkColorSeries, vtkNamedColors
 
 # TODO: Check that these imports are not used -- vtkLookupTable, vtkMinimalStandardRandomSequence, vtkPoints, vtkUnsignedCharArray
-from vtkmodules.vtkCommonComputationalGeometry import vtkParametricTorus
 from vtkmodules.vtkCommonCore import (
     mutable,
     vtkDoubleArray,
     vtkIdList,
     vtkLookupTable,
-    vtkMath,
     vtkPoints,
     vtkUnsignedCharArray,
 )
 from vtkmodules.vtkCommonDataModel import (
-    vtkCellLocator,
     vtkPolyData,
 )
 from vtkmodules.vtkCommonMath import vtkMatrix4x4
 from vtkmodules.vtkCommonTransforms import vtkTransform
 from vtkmodules.vtkFiltersCore import vtkCenterOfMass, vtkGlyph3D, vtkPolyDataNormals
-from vtkmodules.vtkFiltersGeneral import vtkTransformPolyDataFilter
 from vtkmodules.vtkFiltersHybrid import vtkRenderLargeImage
 from vtkmodules.vtkFiltersModeling import vtkBandedPolyDataContourFilter
 from vtkmodules.vtkFiltersSources import (
     vtkArrowSource,
-    vtkDiskSource,
-    vtkLineSource,
-    vtkParametricFunctionSource,
     vtkSphereSource,
 )
 from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
@@ -81,7 +70,7 @@
     vtkPostScriptWriter,
     vtkTIFFWriter,
 )
-from vtkmodules.vtkRenderingAnnotation import vtkAnnotatedCubeActor, vtkAxesActor, vtkScalarBarActor
+from vtkmodules.vtkRenderingAnnotation import vtkAnnotatedCubeActor, vtkAxesActor
 from vtkmodules.vtkRenderingCore import (
     vtkActor,
     vtkPointPicker,
@@ -97,7 +86,6 @@
 import invesalius.constants as const
 import invesalius.data.coordinates as dco
 import invesalius.data.coregistration as dcr
-import invesalius.data.slice_ as sl
 import invesalius.data.styles_3d as styles
 import invesalius.data.transformations as tr
 import invesalius.data.vtk_utils as vtku
@@ -107,11 +95,11 @@
 import invesalius.utils as utils
 from invesalius import inv_paths
 from invesalius.data.actor_factory import ActorFactory
-from invesalius.data.markers.marker import Marker, MarkerType
 from invesalius.data.markers.surface_geometry import SurfaceGeometry
 from invesalius.data.ruler_volume import GenericLeftRulerVolume
 from invesalius.data.visualization.coil_visualizer import CoilVisualizer
 from invesalius.data.visualization.marker_visualizer import MarkerVisualizer
+from invesalius.data.visualization.probe_visualizer import ProbeVisualizer
 from invesalius.data.visualization.vector_field_visualizer import VectorFieldVisualizer
 from invesalius.gui.widgets.canvas_renderer import CanvasRendererCTX
 from invesalius.i18n import tr as _
@@ -326,6 +314,8 @@ def __init__(self, parent):
             vector_field_visualizer=self.vector_field_visualizer,
         )
 
+        self.probe_visualizer = ProbeVisualizer(self.ren, self.interactor)
+
         self.seed_offset = const.SEED_OFFSET
         self.radius_list = vtkIdList()
         self.colors_init = vtkUnsignedCharArray()
@@ -1055,7 +1045,7 @@ def OnUpdateCoilPose(self, m_img, coord):
                 coord[0:3], (self.target_coord[0], -self.target_coord[1], self.target_coord[2])
             )
 
-            formatted_distance = "Distance: {: >5.1f} mm".format(distance_to_target)
+            formatted_distance = f"Distance: {distance_to_target: >5.1f} mm"
 
             if self.distance_text is not None:
                 self.distance_text.SetValue(formatted_distance)
@@ -1235,7 +1225,7 @@ def OnSetTarget(self, marker):
 
         self.coil_visualizer.AddTargetCoil(self.m_target)
 
-        print("Target updated to coordinates {}".format(coord))
+        print(f"Target updated to coordinates {coord}")
 
     def CreateVTKObjectMatrix(self, direction, orientation):
         m_img = dco.coordinates_to_transformation_matrix(
@@ -1504,12 +1494,8 @@ def AddCortexMarkerActor(self, position_orientation, marker_id):
 
             proj = prj.Project()
             timestamp = time.localtime(time.time())
-            stamp_date = "{:0>4d}{:0>2d}{:0>2d}".format(
-                timestamp.tm_year, timestamp.tm_mon, timestamp.tm_mday
-            )
-            stamp_time = "{:0>2d}{:0>2d}{:0>2d}".format(
-                timestamp.tm_hour, timestamp.tm_min, timestamp.tm_sec
-            )
+            stamp_date = f"{timestamp.tm_year:0>4d}{timestamp.tm_mon:0>2d}{timestamp.tm_mday:0>2d}"
+            stamp_time = f"{timestamp.tm_hour:0>2d}{timestamp.tm_min:0>2d}{timestamp.tm_sec:0>2d}"
             sep = "-"
 
             if self.path_meshes is None:
@@ -1567,7 +1553,7 @@ def CreateEfieldSpreadLegend(self):
     def CalculateDistanceMaxEfieldCoGE(self):
         self.distance_efield = distance.euclidean(self.center_gravity_position, self.position_max)
         self.SpreadEfieldFactorTextActor.SetValue(
-            "Spread distance: " + str("{:04.2f}".format(self.distance_efield))
+            "Spread distance: " + str(f"{self.distance_efield:04.2f}")
         )
 
     def EfieldVectors(self):
@@ -1710,7 +1696,7 @@ def find_and_extract_data(self, csv_filename, target_numbers):
 
         matching_rows = []
 
-        with open(csv_filename, "r") as csvfile:
+        with open(csv_filename) as csvfile:
             csv_reader = csv.reader(csvfile)
             for row in csv_reader:
                 # Extract the first three numbers from the current row
@@ -1984,16 +1970,14 @@ def SetEfieldTargetAtCortex(self, position, orientation):
 
     def ShowEfieldAtCortexTarget(self):
         if self.target_at_cortex is not None:
-            import vtk
-
             index = self.efield_mesh.FindPoint(self.target_at_cortex)
             if index in self.Id_list:
                 cell_number = self.Id_list.index(index)
                 self.EfieldAtTargetLegend.SetValue(
-                    "Efield at Target: " + str("{:04.2f}".format(self.e_field_norms[cell_number]))
+                    "Efield at Target: " + str(f"{self.e_field_norms[cell_number]:04.2f}")
                 )
             else:
-                self.EfieldAtTargetLegend.SetValue("Efield at Target: " + str("{:04.2f}".format(0)))
+                self.EfieldAtTargetLegend.SetValue("Efield at Target: " + str(f"{0:04.2f}"))
 
     def CreateEfieldAtTargetLegend(self):
         if self.EfieldAtTargetLegend is not None:
@@ -2143,7 +2127,7 @@ def UpdateEfieldPointLocation(self, m_img, coord, queue_IDs):
             self.e_field_IDs_queue = queue_IDs
             if self.radius_list.GetNumberOfIds() != 0:
                 if np.all(self.old_coord != coord):
-                    self.e_field_IDs_queue.put_nowait((self.radius_list))
+                    self.e_field_IDs_queue.put_nowait(self.radius_list)
                 self.old_coord = np.array([coord])
         except queue.Full:
             pass
@@ -2215,11 +2199,11 @@ def GetEnorm(self, enorm_data, plot_vector):
 
                     proj = prj.Project()
                     timestamp = time.localtime(time.time())
-                    stamp_date = "{:0>4d}{:0>2d}{:0>2d}".format(
-                        timestamp.tm_year, timestamp.tm_mon, timestamp.tm_mday
+                    stamp_date = (
+                        f"{timestamp.tm_year:0>4d}{timestamp.tm_mon:0>2d}{timestamp.tm_mday:0>2d}"
                     )
-                    stamp_time = "{:0>2d}{:0>2d}{:0>2d}".format(
-                        timestamp.tm_hour, timestamp.tm_min, timestamp.tm_sec
+                    stamp_time = (
+                        f"{timestamp.tm_hour:0>2d}{timestamp.tm_min:0>2d}{timestamp.tm_sec:0>2d}"
                     )
                     sep = "-"
 
@@ -2331,8 +2315,6 @@ def SaveEfieldData(self, filename, plot_efield_vectors, marker_id):
     def SavedAllEfieldData(self, filename):
         import csv
 
-        import invesalius.data.imagedata_utils as imagedata_utils
-
         header = [
             "Marker ID",
             "Enorm cell indexes",
diff --git a/invesalius/data/visualization/coil_visualizer.py b/invesalius/data/visualization/coil_visualizer.py
index 3a9504ea0..338a22e50 100644
--- a/invesalius/data/visualization/coil_visualizer.py
+++ b/invesalius/data/visualization/coil_visualizer.py
@@ -36,14 +36,20 @@ def __init__(self, renderer, interactor, actor_factory, vector_field_visualizer)
         # The vector field visualizer is used to show a vector field relative to the coil.
         self.vector_field_visualizer = vector_field_visualizer
 
+        # Number of coils is by default 1
+        self.n_coils = 1
+
         # The actor for showing the actual coil in the volume viewer.
         self.coil_actor = None
+        self.coil_actors = []
 
         # The actor for showing the center of the actual coil in the volume viewer.
         self.coil_center_actor = None
+        self.coil_center_actors = []
 
         # The actor for showing the target coil in the volume viewer.
         self.target_coil_actor = None
+        self.target_coil_actors = []
 
         # The assembly for showing the vector field relative to the coil in the volume viewer.
         self.vector_field_assembly = self.vector_field_visualizer.CreateVectorFieldAssembly()
@@ -73,6 +79,9 @@ def __init__(self, renderer, interactor, actor_factory, vector_field_visualizer)
 
         self.LoadConfig()
 
+        self.AddCoilActor(self.coil_path)
+        self.ShowCoil(False)
+
         self.__bind_events()
 
     def __bind_events(self):
@@ -80,6 +89,7 @@ def __bind_events(self):
         Publisher.subscribe(self.OnNavigationStatus, "Navigation status")
         Publisher.subscribe(self.TrackObject, "Track object")
         Publisher.subscribe(self.ShowCoil, "Show coil in viewer volume")
+        Publisher.subscribe(self.SetCoilCount, "Set coil count")
         Publisher.subscribe(self.ConfigureCoil, "Configure coil")
         Publisher.subscribe(self.UpdateCoilPose, "Update coil pose")
         Publisher.subscribe(self.UpdateVectorField, "Update vector field")
@@ -163,6 +173,9 @@ def SetCoilAtTarget(self, state):
         self.target_coil_actor.GetProperty().SetDiffuseColor(target_coil_color)
         self.coil_center_actor.GetProperty().SetDiffuseColor(target_coil_color)
 
+    def SetCoilCount(self, n_coils):
+        self.n_coils = n_coils
+
     def RemoveCoilActor(self):
         self.renderer.RemoveActor(self.coil_actor)
         self.renderer.RemoveActor(self.coil_center_actor)
@@ -187,16 +200,11 @@ def OnNavigationStatus(self, nav_status, vis_status):
     def TrackObject(self, enabled):
         self.track_object_pressed = enabled
 
-        if self.coil_path is None:
-            return
-
-        # Remove the previous coil actor if it exists.
-        if self.coil_actor is not None:
-            self.RemoveCoilActor()
-
-        # If enabled, add a new coil actor.
+        # Hide coil vector_field if it exists. LUKATODO: what is this vector field?
         if enabled:
-            self.AddCoilActor(self.coil_path)
+            self.vector_field_assembly.SetVisibility(1)
+        else:
+            self.vector_field_assembly.SetVisibility(0)
 
     # Called when 'show coil' button is pressed in the user interface or in code.
     def ShowCoil(self, state):
diff --git a/invesalius/data/visualization/probe_visualizer.py b/invesalius/data/visualization/probe_visualizer.py
new file mode 100644
index 000000000..3627bf7b3
--- /dev/null
+++ b/invesalius/data/visualization/probe_visualizer.py
@@ -0,0 +1,161 @@
+import os
+
+import vtk
+
+import invesalius.data.vtk_utils as vtku
+from invesalius import inv_paths
+from invesalius.pubsub import pub as Publisher
+
+
+class ProbeVisualizer:
+    """
+    A class for visualizing probe in the volume viewer.
+    """
+
+    def __init__(self, renderer, interactor):
+        self.renderer = renderer
+        self.interactor = interactor
+
+        self.probe_actor = None
+        self.probe_path = os.path.join(inv_paths.OBJ_DIR, "stylus.stl")
+        self.show_probe = False
+        self.is_navigating = False
+
+        self.AddProbeActor(self.probe_path)
+        self.__bind_events()
+
+    def __bind_events(self):
+        Publisher.subscribe(self.ShowProbe, "Show probe in viewer volume")
+        Publisher.subscribe(self.UpdateProbePose, "Update probe pose")
+
+    def ShowProbe(self, state):
+        self.show_probe = state
+
+        if self.probe_actor:
+            self.probe_actor.SetVisibility(self.show_probe)
+            self.interactor.Render()
+
+    def AddProbeActor(self, probe_path):
+        vtk_colors = vtk.vtkNamedColors()
+        obj_polydata = vtku.CreateObjectPolyData(probe_path)
+
+        # This rotation is done so that the probe visibility icon in viewer_volume doesn"t need to be rotated
+        transform = vtk.vtkTransform()
+        transform.RotateZ(-150)
+
+        transform_filt = vtk.vtkTransformPolyDataFilter()
+        transform_filt.SetTransform(transform)
+        transform_filt.SetInputData(obj_polydata)
+        transform_filt.Update()
+
+        obj_mapper = vtk.vtkPolyDataMapper()
+        obj_mapper.SetInputData(transform_filt.GetOutput())
+        obj_mapper.ScalarVisibilityOff()
+
+        probe_actor = vtk.vtkActor()
+        probe_actor.SetMapper(obj_mapper)
+        probe_actor.GetProperty().SetAmbientColor(vtk_colors.GetColor3d("GhostWhite"))
+        probe_actor.GetProperty().SetSpecular(30)
+        probe_actor.GetProperty().SetSpecularPower(80)
+        probe_actor.GetProperty().SetOpacity(0.4)
+        probe_actor.SetVisibility(0)
+
+        self.probe_actor = probe_actor
+        self.renderer.AddActor(self.probe_actor)
+
+    def RemoveProbeActor(self):
+        if self.probe_actor is not None:
+            self.renderer.RemoveActor(self.probe_actor)
+            self.probe_actor = None
+
+    def UpdateProbePose(self, m_img, coord):
+        m_img_flip = m_img.copy()
+        m_img_flip[1, -1] = -m_img_flip[1, -1]
+
+        m_img_vtk = vtku.numpy_to_vtkMatrix4x4(m_img_flip)
+
+        self.probe_actor.SetUserMatrix(m_img_vtk)
+
+
+"""
+Code used to shrink and align old stylus.stl file. 
+May be useful for aligning other stl files into VTK in the future.
+This computes the main axes of the object and position of the tip
+then rotates and translates to align with vtk axes 
+
+    import numpy as np
+
+    import invesalius.data.coordinates as dco
+    import invesalius.constants as const
+    import invesalius.data.polydata_utils as pu
+
+    obj_polydata = vtku.CreateObjectPolyData(probe_path) # STL file to be aligned
+    #Shrink object
+    scale_factor = 0.1
+    shrink_transform = vtk.vtkTransform()
+    shrink_transform.Scale(scale_factor, scale_factor, scale_factor)
+
+    shrink_transform_filt = vtk.vtkTransformPolyDataFilter()
+    shrink_transform_filt.SetTransform(shrink_transform)
+    shrink_transform_filt.SetInputData(obj_polydata)
+    shrink_transform_filt.Update()
+
+    obj_polydata = shrink_transform_filt.GetOutput() #replace original with shrunk
+
+    # Rotate to align obj_axes and translate tip to origin
+    obj_axes, tip = self.__calculate_pca(obj_polydata)
+
+    obj_to_vtk = np.linalg.inv(obj_axes) 
+    rotation = vtk.vtkMatrix4x4() 
+    rotation.SetElement(3, 3, 1) #affine transformation: last row is [0 0 0 1]
+    for i in range(3): # copy obj_to_vtk to vtk_matrix
+        for j in range(3):
+            rotation.SetElement(i, j, obj_to_vtk[i, j])
+            rotation.SetElement(3, j, 0) #last column
+        rotation.SetElement(i, 3, 0) # last row
+
+    transform = vtk.vtkTransform()
+    transform.Translate(-obj_to_vtk@tip) #change to vtk basis before translating
+    transform.Concatenate(rotation)
+
+    transform_filt = vtk.vtkTransformPolyDataFilter()
+    transform_filt.SetTransform(transform)
+    transform_filt.SetInputData(obj_polydata)
+    transform_filt.Update()
+
+    # Write the aligned object into a new STL file
+    output_path = os.path.join(inv_paths.OBJ_DIR, "stylus2.stl")
+    stl_writer = vtk.vtkSTLWriter()
+    stl_writer.SetFileTypeToBinary()
+    stl_writer.SetFileName(output_path)
+    stl_writer.SetInputData(transform_filt.GetOutput())
+    stl_writer.Write()
+
+    # Helper function
+    def __calculate_pca(self, polydata): 
+        #Calculate 3 principal axes of object from vtk polydata with principal component analysis.
+        #Returns the object axes in the desired xyz-order in a 3x3 matrix:
+        #    x-column: Secondary axis (lateral axis of object)
+        #    y-column: Principal axis (longitudinal axis of object)
+        #    z-column: Tertiary axis (normal to object)
+        #Google aircraft axes for analogue
+
+        points = np.array(polydata.GetPoints().GetData())
+
+        # Center the points
+        centroid = np.mean(points, axis=0)
+        centered_points = points - centroid
+
+        # Perform PCA
+        covariance_matrix = np.cov(centered_points, rowvar=False)
+        eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)
+        #eigh returns vectors in ascending order (ie. principal axis is the last one)
+        
+        # Permute, so that y-column is the principal axis and x-column secondary axis
+        eigenvectors[:, [0, 1, 2]] = eigenvectors[:, [1, 2, 0]]
+        
+        # the tip is the point furthest along principal axis 
+        # for another STL object tip could be in opposite direction so would use argmin instead
+        tip = points[np.argmax(np.dot(centered_points, eigenvectors[:,1]))]
+        return eigenvectors, tip
+"""
diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py
index b84201484..8cdf788f6 100644
--- a/invesalius/gui/dialogs.py
+++ b/invesalius/gui/dialogs.py
@@ -4461,9 +4461,7 @@ def RemoveSinglePointActor(self):
 
     def GetCurrentCoord(self):
         coord_raw, marker_visibilities = self.tracker.TrackerCoordinates.GetCoordinates()
-        coord, _ = dcr.corregistrate_dynamic(
-            (self.m_change, 0), coord_raw, const.DEFAULT_REF_MODE, [None, None]
-        )
+        coord, _ = dcr.corregistrate_probe(self.m_change, None, coord_raw, const.DEFAULT_REF_MODE)
         return coord[:3], marker_visibilities
 
     def AddMarker(self, size, colour, coord):
diff --git a/invesalius/gui/preferences.py b/invesalius/gui/preferences.py
index 18fc9c057..f422951f6 100644
--- a/invesalius/gui/preferences.py
+++ b/invesalius/gui/preferences.py
@@ -947,10 +947,7 @@ def OnChooseTracker(self, evt, ctrl):
             self.ShowParent()
 
     def OnChooseReferenceMode(self, evt, ctrl):
-        # Probably need to refactor object registration as a whole to use the
-        # OnChooseReferenceMode function which was used earlier. It can be found in
-        # the deprecated code in ObjectRegistrationPanel in task_navigator.py.
-        pass
+        Navigation().SetReferenceMode(evt.GetSelection())
 
     def HideParent(self):  # hide preferences dialog box
         self.GetGrandParent().Hide()
diff --git a/invesalius/gui/task_navigator.py b/invesalius/gui/task_navigator.py
index 8279db5c4..1d8d50dda 100644
--- a/invesalius/gui/task_navigator.py
+++ b/invesalius/gui/task_navigator.py
@@ -327,6 +327,7 @@ def __init__(self, parent, nav_hub):
         book.AddPage(ImagePage(book, nav_hub), _("Image"))
         book.AddPage(TrackerPage(book, nav_hub), _("Patient"))
         book.AddPage(RefinePage(book, nav_hub), _("Refine"))
+        book.AddPage(StylusPage(book, nav_hub), _("Stylus"))
         book.AddPage(StimulatorPage(book, nav_hub), _("TMS Coil"))
 
         book.SetSelection(0)
@@ -342,6 +343,7 @@ def __init__(self, parent, nav_hub):
     def __bind_events(self):
         Publisher.subscribe(self._FoldTracker, "Move to tracker page")
         Publisher.subscribe(self._FoldRefine, "Move to refine page")
+        Publisher.subscribe(self._FoldStylus, "Move to stylus page")
         Publisher.subscribe(self._FoldStimulator, "Move to stimulator page")
         Publisher.subscribe(self._FoldImage, "Move to image page")
 
@@ -363,7 +365,7 @@ def OnPageChanged(self, evt):
             Publisher.sendMessage("Update UI for refine tab")
 
         # new page validations
-        if (old_page == 1) and (new_page == 2 or new_page == 3):
+        if (old_page == 1) and (new_page > 1):
             # Do not allow user to move to other (forward) tabs if tracker fiducials not done.
             if self.image.AreImageFiducialsSet() and not self.tracker.AreTrackerFiducialsSet():
                 self.book.SetSelection(1)
@@ -380,9 +382,12 @@ def _FoldTracker(self):
     def _FoldRefine(self):
         self.book.SetSelection(2)
 
-    def _FoldStimulator(self):
+    def _FoldStylus(self):
         self.book.SetSelection(3)
 
+    def _FoldStimulator(self):
+        self.book.SetSelection(4)
+
 
 class ImagePage(wx.Panel):
     def __init__(self, parent, nav_hub):
@@ -508,9 +513,7 @@ def LoadImageFiducials(self, label, position):
     def GetFiducialByAttribute(self, fiducials, attribute_name, attribute_value):
         found = [fiducial for fiducial in fiducials if fiducial[attribute_name] == attribute_value]
 
-        assert len(found) != 0, "No fiducial found for which {} = {}".format(
-            attribute_name, attribute_value
-        )
+        assert len(found) != 0, f"No fiducial found for which {attribute_name} = {attribute_value}"
         return found[0]
 
     def SetImageFiducial(self, fiducial_name, position):
@@ -764,9 +767,7 @@ def StopRegistration(self):
     def GetFiducialByAttribute(self, fiducials, attribute_name, attribute_value):
         found = [fiducial for fiducial in fiducials if fiducial[attribute_name] == attribute_value]
 
-        assert len(found) != 0, "No fiducial found for which {} = {}".format(
-            attribute_name, attribute_value
-        )
+        assert len(found) != 0, f"No fiducial found for which {attribute_name} = {attribute_value}"
         return found[0]
 
     def OnSetTrackerFiducial(self, fiducial_name):
@@ -997,7 +998,7 @@ def OnBack(self, evt):
         Publisher.sendMessage("Move to tracker page")
 
     def OnNext(self, evt):
-        Publisher.sendMessage("Move to stimulator page")
+        Publisher.sendMessage("Move to stylus page")
 
     def OnRefine(self, evt):
         self.icp.RegisterICP(self.navigation, self.tracker)
@@ -1005,6 +1006,102 @@ def OnRefine(self, evt):
             self.OnUpdateUI()
 
 
+class StylusPage(wx.Panel):
+    def __init__(self, parent, nav_hub):
+        wx.Panel.__init__(self, parent)
+        self.navigation = nav_hub.navigation
+        self.tracker = nav_hub.tracker
+
+        border = wx.FlexGridSizer(1, 3, 5)
+        self.border = border
+
+        self.done = False
+
+        lbl = wx.StaticText(self, -1, _("Calibrate stylus with head"))
+        lbl.SetFont(wx.Font(9, wx.DEFAULT, wx.NORMAL, wx.BOLD))
+        self.lbl = lbl
+        self.help_img = wx.Image(os.path.join(inv_paths.ICON_DIR, "align.png"), wx.BITMAP_TYPE_ANY)
+
+        # first show help in grayscale. when record successful: make it green to show success
+        self.help = wx.GenericStaticBitmap(
+            self,
+            -1,
+            self.help_img.ConvertToGreyscale(),
+            (10, 5),
+            (self.help_img.GetWidth(), self.help_img.GetHeight()),
+        )
+
+        lbl_rec = wx.StaticText(self, -1, _("Point stylus up relative to head, like so:"))
+        btn_rec = wx.Button(self, -1, _("Record"))
+        btn_rec.SetToolTip("Record stylus orientation relative to head")
+        btn_rec.Bind(wx.EVT_BUTTON, self.onRecord)
+
+        self.btn_rec = btn_rec
+
+        self.border.AddMany(
+            [
+                (lbl, 1, wx.EXPAND | wx.ALL | wx.ALIGN_CENTER_VERTICAL, 2),
+                (lbl_rec, 1, wx.EXPAND | wx.ALL | wx.ALIGN_CENTER_VERTICAL, 2),
+                (self.help, 0, wx.EXPAND | wx.ALL | wx.ALIGN_CENTER_VERTICAL, 1),
+                (btn_rec, 0, wx.EXPAND | wx.ALL | wx.ALIGN_LEFT, 1),
+            ]
+        )
+
+        back_button = wx.Button(self, label="Back")
+        back_button.Bind(wx.EVT_BUTTON, partial(self.OnBack))
+        next_button = wx.Button(self, label="Next")
+        next_button.Bind(wx.EVT_BUTTON, partial(self.OnNext))
+
+        bottom_sizer = wx.BoxSizer(wx.HORIZONTAL)
+        bottom_sizer.Add(back_button)
+        bottom_sizer.Add(next_button)
+
+        main_sizer = wx.BoxSizer(wx.VERTICAL)
+        main_sizer.AddMany(
+            [
+                (border, 0, wx.ALIGN_CENTER_HORIZONTAL | wx.ALL, 1),
+                (bottom_sizer, 0, wx.ALIGN_CENTER | wx.CENTER | wx.TOP, 1),
+            ]
+        )
+
+        self.SetSizerAndFit(main_sizer)
+        self.Layout()
+        self.__bind_events()
+
+    def __bind_events(self):
+        pass  # Publisher.subscribe(self.OnCloseProject, 'Remove object data')
+
+    def onRecord(self, evt):
+        marker_visibilities, __, coord_raw = self.tracker.GetTrackerCoordinates(
+            ref_mode_id=0, n_samples=1
+        )
+
+        if marker_visibilities[0] and marker_visibilities[1]:  # if probe and head are visible
+            if self.navigation.SetStylusOrientation(coord_raw) and not self.done:
+                # if successfully created r_stylus in navigation for the first time: make the illustration green to show success
+                self.done = True
+                self.help.Destroy()  # show a colored (green) bitmap as opposed to grayscale
+                self.help = wx.GenericStaticBitmap(
+                    self,
+                    -1,
+                    self.help_img.ConvertToBitmap(),
+                    (10, 5),
+                    (self.help_img.GetWidth(), self.help_img.GetHeight()),
+                )
+                self.border.Insert(
+                    2, self.help, 0, wx.EXPAND | wx.ALL | wx.ALIGN_CENTER_VERTICAL, 1
+                )
+                self.Layout()
+        else:
+            wx.MessageBox(_("Probe or head not visible to tracker!"), _("InVesalius 3"))
+
+    def OnBack(self, evt):
+        Publisher.sendMessage("Move to refine page")
+
+    def OnNext(self, evt):
+        Publisher.sendMessage("Move to stimulator page")
+
+
 class StimulatorPage(wx.Panel):
     def __init__(self, parent, nav_hub):
         wx.Panel.__init__(self, parent)
@@ -1207,8 +1304,6 @@ def __init__(self, parent, nav_hub):
         track_object_button.SetBackgroundColour(GREY_COLOR)
         track_object_button.SetBitmap(BMP_TRACK)
         track_object_button.SetValue(False)
-        if not self.track_obj:
-            track_object_button.Enable(False)
         track_object_button.SetToolTip(tooltip)
         track_object_button.Bind(
             wx.EVT_TOGGLEBUTTON, partial(self.OnTrackObjectButton, ctrl=track_object_button)
@@ -1243,10 +1338,26 @@ def __init__(self, parent, nav_hub):
         show_coil_button.SetBitmap(BMP_SHOW_COIL)
         show_coil_button.SetToolTip(tooltip)
         show_coil_button.SetValue(False)
-        show_coil_button.Enable(False)
+        show_coil_button.Enable(True)
         show_coil_button.Bind(wx.EVT_TOGGLEBUTTON, self.OnShowCoil)
         self.show_coil_button = show_coil_button
 
+        # Toggle button for showing probe during navigation
+        tooltip = _("Show probe")
+        BMP_SHOW_PROBE = wx.Bitmap(
+            str(inv_paths.ICON_DIR.joinpath("stylus.png")), wx.BITMAP_TYPE_PNG
+        )
+        show_probe_button = wx.ToggleButton(
+            self, -1, "", style=pbtn.PB_STYLE_SQUARE, size=ICON_SIZE
+        )
+        show_probe_button.SetBackgroundColour(GREY_COLOR)
+        show_probe_button.SetBitmap(BMP_SHOW_PROBE)
+        show_probe_button.SetToolTip(tooltip)
+        show_probe_button.Enable(True)
+        self.UpdateToggleButton(show_probe_button, False)  # the probe is hidden at start
+        show_probe_button.Bind(wx.EVT_TOGGLEBUTTON, self.OnShowProbe)
+        self.show_probe_button = show_probe_button
+
         # Toggle Button to use serial port to trigger pulse signal and create markers
         tooltip = _("Enable serial port communication to trigger pulse and create markers")
         BMP_PORT = wx.Bitmap(str(inv_paths.ICON_DIR.joinpath("wave.png")), wx.BITMAP_TYPE_PNG)
@@ -1347,6 +1458,7 @@ def __init__(self, parent, nav_hub):
                 (efield_checkbox),
                 (lock_to_target_button),
                 (show_coil_button),
+                (show_probe_button),
             ]
         )
 
@@ -1388,6 +1500,8 @@ def __bind_events(self):
         Publisher.subscribe(self.UpdateTractsVisualization, "Update tracts visualization")
 
         # Externally press/unpress and enable/disable buttons.
+        Publisher.subscribe(self.PressShowProbeButton, "Press show-probe button")
+
         Publisher.subscribe(self.PressShowCoilButton, "Press show-coil button")
         Publisher.subscribe(self.EnableShowCoilButton, "Enable show-coil button")
 
@@ -1551,12 +1665,11 @@ def OnCheckStatus(self, nav_status, vis_status):
             self.EnableToggleButton(self.checkbox_serial_port, 1)
             self.UpdateToggleButton(self.checkbox_serial_port)
 
-        # Enable/Disable track-object checkbox if navigation is off/on and object registration is valid.
+        # Enable/Disable track-object checkbox if object registration is valid.
         obj_registration = self.navigation.GetObjectRegistration()
-        enable_track_object = (
-            obj_registration is not None and obj_registration[0] is not None and not nav_status
-        )
+        enable_track_object = obj_registration is not None and obj_registration[0] is not None
         self.EnableTrackObjectButton(enable_track_object)
+        self.EnableShowCoilButton(enable_track_object)
 
     # Robot
     def OnRobotStatus(self, data):
@@ -1653,11 +1766,9 @@ def OnTrackObjectButton(self, evt=None, ctrl=None):
         if not pressed:
             Publisher.sendMessage("Press target mode button", pressed=pressed)
 
-        # Disable or enable 'Show coil' button, based on if 'Track object' button is pressed.
-        Publisher.sendMessage("Enable show-coil button", enabled=pressed)
-
-        # Also, automatically press or unpress 'Show coil' button.
+        # Automatically press or unpress 'Show coil' and 'Show probe' button.
         Publisher.sendMessage("Press show-coil button", pressed=pressed)
+        Publisher.sendMessage("Press show-probe button", pressed=(not pressed))
 
         self.SaveConfig()
 
@@ -1681,6 +1792,16 @@ def OnShowCoil(self, evt=None):
         pressed = self.show_coil_button.GetValue()
         Publisher.sendMessage("Show coil in viewer volume", state=pressed)
 
+    # 'Show probe' button
+    def PressShowProbeButton(self, pressed=False):
+        self.UpdateToggleButton(self.show_probe_button, pressed)
+        self.OnShowProbe()
+
+    def OnShowProbe(self, evt=None):
+        self.UpdateToggleButton(self.show_probe_button)
+        pressed = self.show_probe_button.GetValue()
+        Publisher.sendMessage("Show probe in viewer volume", state=pressed)
+
     # 'Serial Port Com'
     def OnEnableSerialPort(self, evt, ctrl):
         self.UpdateToggleButton(ctrl)
@@ -2852,7 +2973,7 @@ def ParseValue(self, value):
 
     def GetMarkersFromFile(self, filename, overwrite_image_fiducials):
         try:
-            with open(filename, "r") as file:
+            with open(filename) as file:
                 magick_line = file.readline()
                 assert magick_line.startswith(const.MARKER_FILE_MAGICK_STRING)
                 version = int(magick_line.split("_")[-1])
@@ -2926,12 +3047,8 @@ def OnShowHideAllMarkers(self, evt, ctrl):
     def OnSaveMarkers(self, evt):
         prj_data = prj.Project()
         timestamp = time.localtime(time.time())
-        stamp_date = "{:0>4d}{:0>2d}{:0>2d}".format(
-            timestamp.tm_year, timestamp.tm_mon, timestamp.tm_mday
-        )
-        stamp_time = "{:0>2d}{:0>2d}{:0>2d}".format(
-            timestamp.tm_hour, timestamp.tm_min, timestamp.tm_sec
-        )
+        stamp_date = f"{timestamp.tm_year:0>4d}{timestamp.tm_mon:0>2d}{timestamp.tm_mday:0>2d}"
+        stamp_time = f"{timestamp.tm_hour:0>2d}{timestamp.tm_min:0>2d}{timestamp.tm_sec:0>2d}"
         sep = "-"
         parts = [stamp_date, stamp_time, prj_data.name, "markers"]
         default_filename = sep.join(parts) + ".mkss"
diff --git a/invesalius/navigation/navigation.py b/invesalius/navigation/navigation.py
index 54d8e87ee..ca08e51f2 100644
--- a/invesalius/navigation/navigation.py
+++ b/invesalius/navigation/navigation.py
@@ -35,7 +35,6 @@
 import invesalius.data.vtk_utils as vtk_utils
 import invesalius.project as prj
 import invesalius.session as ses
-import invesalius.utils as utils
 from invesalius.data.markers.marker import MarkerType
 from invesalius.navigation.image import Image
 from invesalius.navigation.iterativeclosestpoint import IterativeClosestPoint
@@ -133,15 +132,28 @@ def __init__(self, vis_queues, vis_components, event, sle, neuronavigation_api):
         self.sle = sle
         self.event = event
         self.neuronavigation_api = neuronavigation_api
+        self.navigation = Navigation()
 
     def run(self):
         while not self.event.is_set():
             got_coords = False
-            object_visible_flag = False
             try:
-                coord, marker_visibilities, m_img, view_obj = self.coord_queue.get_nowait()
+                coords, marker_visibilities, m_imgs = self.coord_queue.get_nowait()
                 got_coords = True
-                object_visible_flag = marker_visibilities[2]
+
+                probe_visible = marker_visibilities[0]
+                coil_visible = marker_visibilities[2]
+
+                # automatically track either coil or stylus if only one of them is visible, otherise use navigation.track_obj
+                track_coil = (
+                    (coil_visible or not probe_visible)
+                    if (coil_visible ^ probe_visible)
+                    else self.navigation.track_obj
+                )
+
+                # choose which object to track in slices and viewer_volume pointer
+                coord = coords[track_coil]  # main-coil if track_coil else stylus
+                m_img = m_imgs[track_coil]
 
                 # use of CallAfter is mandatory otherwise crashes the wx interface
                 if self.view_tracts:
@@ -178,15 +190,32 @@ def run(self):
                 # new marker is created, it is created in the current position of the object.
                 wx.CallAfter(Publisher.sendMessage, "Set cross focal point", position=coord)
 
-                if self.e_field_loaded and object_visible_flag:
+                wx.CallAfter(
+                    Publisher.sendMessage,
+                    "Update volume viewer pointer",
+                    position=[coord[0], -coord[1], coord[2]],
+                )
+
+                if coil_visible:
+                    wx.CallAfter(
+                        Publisher.sendMessage, "Update coil pose", m_img=m_imgs[1], coord=coords[1]
+                    )
                     wx.CallAfter(
                         Publisher.sendMessage,
-                        "Update point location for e-field calculation",
-                        m_img=m_img,
-                        coord=coord,
-                        queue_IDs=self.e_field_IDs_queue,
+                        "Update object arrow matrix",
+                        m_img=m_imgs[1],
+                        coord=coords[1],
+                        flag=self.peel_loaded,
                     )
-                    if not self.e_field_norms_queue.empty():
+
+                    if self.e_field_loaded:
+                        wx.CallAfter(
+                            Publisher.sendMessage,
+                            "Update point location for e-field calculation",
+                            m_img=m_imgs[1],
+                            coord=coords[1],
+                            queue_IDs=self.e_field_IDs_queue,
+                        )
                         try:
                             enorm_data = self.e_field_norms_queue.get_nowait()
                             wx.CallAfter(
@@ -195,26 +224,16 @@ def run(self):
                                 enorm_data=enorm_data,
                                 plot_vector=self.plot_efield_vectors,
                             )
-                        finally:
+                        except queue.Empty:
+                            pass
+                        else:
                             self.e_field_norms_queue.task_done()
 
-                if view_obj:
-                    wx.CallAfter(
-                        Publisher.sendMessage, "Update coil pose", m_img=m_img, coord=coord
-                    )
+                if probe_visible:
                     wx.CallAfter(
-                        Publisher.sendMessage,
-                        "Update object arrow matrix",
-                        m_img=m_img,
-                        coord=coord,
-                        flag=self.peel_loaded,
-                    )
-                else:
-                    wx.CallAfter(
-                        Publisher.sendMessage,
-                        "Update volume viewer pointer",
-                        position=[coord[0], -coord[1], coord[2]],
+                        Publisher.sendMessage, "Update probe pose", m_img=m_imgs[0], coord=coords[0]
                     )
+
                 # Render the volume viewer and the slice viewers.
                 wx.CallAfter(Publisher.sendMessage, "Render volume viewer")
                 wx.CallAfter(Publisher.sendMessage, "Update slice viewer")
@@ -235,9 +254,12 @@ def __init__(self, pedal_connector, neuronavigation_api):
 
         self.correg = None
         self.target = None
+        self.n_coils = 1
+        self.obj_registrations = []
         self.object_registration = None
         self.track_obj = False
         self.m_change = None
+        self.r_stylus = None
         self.obj_data = None
         self.all_fiducials = np.zeros((6, 6))
         self.event = threading.Event()
@@ -322,19 +344,21 @@ def LoadConfig(self):
         session = ses.Session()
         state = session.GetConfig("navigation")
 
-        if state is None:
-            return
+        if state is not None:
+            object_fiducials = np.array(state["object_fiducials"])
+            object_orientations = np.array(state["object_orientations"])
+            object_reference_mode = state["object_reference_mode"]
+            object_name = state["object_name"].encode(const.FS_ENCODE)
+            self.object_registration = (
+                object_fiducials,
+                object_orientations,
+                object_reference_mode,
+                object_name,
+            )
 
-        object_fiducials = np.array(state["object_fiducials"])
-        object_orientations = np.array(state["object_orientations"])
-        object_reference_mode = state["object_reference_mode"]
-        object_name = state["object_name"].encode(const.FS_ENCODE)
-        self.object_registration = (
-            object_fiducials,
-            object_orientations,
-            object_reference_mode,
-            object_name,
-        )
+            # Try to load stylus orientation data
+            if "r_stylus" in state:
+                self.r_stylus = np.array(state["r_stylus"])
 
     def CoilAtTarget(self, state):
         self.coil_at_target = state
@@ -403,6 +427,35 @@ def EstimateTrackerToInVTransformationMatrix(self, tracker, image):
             self.all_fiducials[3:, :].T, self.all_fiducials[:3, :].T, shear=False, scale=False
         )
 
+    def SetStylusOrientation(self, coord_raw):
+        if self.m_change is not None:
+            m_probe = dcr.compute_marker_transformation(coord_raw, 0)
+
+            # transform probe to reference system if dynamic ref_mode
+            if self.ref_mode_id:
+                up_trk = dcr.object_to_reference(coord_raw, m_probe)[:3, :3]
+            else:
+                up_trk = m_probe[:3, :3]
+
+            # up_trk: orientation 'stylus pointing up along head' in tracker space
+            # up_vtk: orientation 'stylus pointing up along head' in vtk space
+            up_vtk = tr.euler_matrix(*np.radians([90.0, 0.0, 0.0]), axes="rxyz")[:3, :3]
+
+            # Rotate 90 degrees around z-axis from trk system where stylus points in x-axis
+            # to vtk-system where stylus points in y-axis
+            R = tr.euler_matrix(*np.radians([0, 0, -90]), axes="rxyz")[:3, :3]
+
+            # Rotation from tracker to VTK coordinate system
+            self.r_stylus = up_vtk @ np.linalg.inv(R @ up_trk @ np.linalg.inv(R))
+            # Save r_stylus to config file
+            session = ses.Session()
+            if (nav_config := session.GetConfig("navigation")) is not None:
+                nav_config["r_stylus"] = self.r_stylus.tolist()
+                session.SetConfig("navigation", nav_config)
+            return True
+        else:
+            return False
+
     def StartNavigation(self, tracker, icp):
         # initialize jobs list
         jobs_list = []
@@ -427,71 +480,53 @@ def StartNavigation(self, tracker, icp):
         ]
 
         Publisher.sendMessage("Navigation status", nav_status=True, vis_status=vis_components)
-        errors = False
-
-        if self.track_obj:
-            # if object tracking is selected
-            if self.object_registration is None:
-                # 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
-                # object_registration[0] is object 3x3 fiducial matrix and object_registration[1] is 3x3 orientation matrix
-                obj_fiducials, obj_orients, obj_ref_mode, obj_name = self.object_registration
 
-                coreg_data = [self.m_change, obj_ref_mode]
+        # LUKATODO: object_registration --> coil_registrations
+        if self.object_registration is None:
+            # check if object registration was performed
+            wx.MessageBox(_("Perform coil registration before navigation."), _("InVesalius 3"))
+        else:
+            # if object registration was correctly performed continue with navigation
+            # object_registration[0] is object 3x3 fiducial matrix and object_registration[1] is 3x3 orientation matrix
+            obj_fiducials, obj_orients, obj_ref_mode, obj_name = self.object_registration
+
+            coreg_data = [self.r_stylus, self.m_change, obj_ref_mode]
 
-                if self.ref_mode_id:
-                    coord_raw, marker_visibilities = tracker.TrackerCoordinates.GetCoordinates()
-                else:
-                    coord_raw = np.array([None])
+            if self.ref_mode_id:
+                coord_raw, marker_visibilities = tracker.TrackerCoordinates.GetCoordinates()
+            else:
+                coord_raw = np.array([None])
 
-                self.obj_data = db.object_registration(
-                    obj_fiducials, obj_orients, coord_raw, self.m_change
-                )
-                coreg_data.extend(self.obj_data)
-
-                queues = [
-                    self.coord_queue,
-                    self.coord_tracts_queue,
-                    self.icp_queue,
-                    self.object_at_target_queue,
-                    self.efield_queue,
-                ]
-                jobs_list.append(
-                    dcr.CoordinateCorregistrate(
-                        self.ref_mode_id,
-                        tracker,
-                        coreg_data,
-                        self.view_tracts,
-                        queues,
-                        self.event,
-                        self.sleep_nav,
-                        tracker.tracker_id,
-                        self.target,
-                        icp,
-                        self.e_field_loaded,
-                    )
-                )
-        else:
-            coreg_data = (self.m_change, 0)
-            queues = [self.coord_queue, self.coord_tracts_queue, self.icp_queue, self.efield_queue]
+            # LUKATODO: NOTE: coil does not need to be visible here, we only need the head coord for obj_data:
+            self.obj_data = db.object_registration(
+                obj_fiducials, obj_orients, coord_raw, self.m_change
+            )
+            coreg_data.extend(self.obj_data)
+
+            queues = [
+                self.coord_queue,
+                self.coord_tracts_queue,
+                self.icp_queue,
+                self.object_at_target_queue,
+                self.efield_queue,
+            ]
             jobs_list.append(
-                dcr.CoordinateCorregistrateNoObject(
+                dcr.CoordinateCorregistrate(
                     self.ref_mode_id,
                     tracker,
+                    self.n_coils,
                     coreg_data,
                     self.view_tracts,
                     queues,
                     self.event,
                     self.sleep_nav,
+                    tracker.tracker_id,
+                    self.target,
                     icp,
                     self.e_field_loaded,
                 )
             )
 
-        if not errors:
             # TODO: Test the serial port thread
             if self.serial_port_in_use:
                 self.serial_port_connection = spc.SerialPortConnection(
diff --git a/navigation/ndi_files/Markers/NBSref.rom b/navigation/ndi_files/Markers/NBSref.rom
index 426497352..820d7b8fb 100644
Binary files a/navigation/ndi_files/Markers/NBSref.rom and b/navigation/ndi_files/Markers/NBSref.rom differ
diff --git a/navigation/ndi_files/Markers/NT-115.rom b/navigation/ndi_files/Markers/NT-115.rom
index dd67e7678..a5868cbad 100644
Binary files a/navigation/ndi_files/Markers/NT-115.rom and b/navigation/ndi_files/Markers/NT-115.rom differ
diff --git a/navigation/objects/stylus.stl b/navigation/objects/stylus.stl
index bf8058efb..e94aa6b32 100644
Binary files a/navigation/objects/stylus.stl and b/navigation/objects/stylus.stl differ