coordinates.py 10.1 KB
#--------------------------------------------------------------------------
# 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.
#--------------------------------------------------------------------------

from math import sin, cos
import numpy as np

from time import sleep
from random import uniform


def GetCoordinates(trck_init, trck_id, ref_mode):

    """
    Read coordinates from spatial tracking devices using

    :param trck_init: Initialization variable of tracking device and connection type. See tracker.py.
    :param trck_id: ID of tracking device.
    :param ref_mode: Single or dynamic reference mode of tracking.
    :return: array of six coordinates (x, y, z, alpha, beta, gamma)
    """

    coord = None
    if trck_id:
        getcoord = {1: ClaronCoord,
                    2: PolhemusCoord,
                    3: PolhemusCoord,
                    4: PolhemusCoord,
                    5: DebugCoord}
        coord = getcoord[trck_id](trck_init, trck_id, ref_mode)
    else:
        print "Select Tracker"

    return coord


def ClaronCoord(trck_init, trck_id, ref_mode):
    trck = trck_init[0]
    scale = np.array([1.0, 1.0, -1.0])
    coord = None
    k = 0
    # TODO: try to replace while and use some Claron internal computation
    if ref_mode:
        while k < 20:
            try:
                trck.Run()
                probe = np.array([trck.PositionTooltipX1 * scale[0], trck.PositionTooltipY1 * scale[1],
                                  trck.PositionTooltipZ1 * scale[2], trck.AngleX1, trck.AngleY1, trck.AngleZ1])
                reference = np.array([trck.PositionTooltipX2 * scale[0], trck.PositionTooltipY2 * scale[1],
                                      trck.PositionTooltipZ2 * scale[2], trck.AngleX2, trck.AngleY2, trck.AngleZ2])
                k = 30
            except AttributeError:
                k += 1
                print "wait, collecting coordinates ..."
        if k == 30:
            coord = dynamic_reference(probe, reference)
    else:
        while k < 20:
            try:
                trck.Run()
                coord = np.array([trck.PositionTooltipX1 * scale[0], trck.PositionTooltipY1 * scale[1],
                                  trck.PositionTooltipZ1 * scale[2], trck.AngleX1, trck.AngleY1, trck.AngleZ1])
                k = 30
            except AttributeError:
                k += 1
                print "wait, collecting coordinates ..."

    return coord


def PolhemusCoord(trck, trck_id, ref_mode):
    coord = None

    if trck[1] == 'serial':
        coord = PolhemusSerialCoord(trck[0], trck_id, ref_mode)

    elif trck[1] == 'usb':
        coord = PolhemusUSBCoord(trck[0], trck_id, ref_mode)

    elif trck[1] == 'wrapper':
        coord = PolhemusWrapperCoord(trck[0], trck_id, ref_mode)

    return coord


def PolhemusWrapperCoord(trck, trck_id, ref_mode):
    if trck_id == 2:
        scale = 10. * np.array([1., 1.0, -1.0])
    else:
        scale = 25.4 * np.array([1., 1.0, -1.0])
    coord = None

    if ref_mode:
        trck.Run()
        probe = np.array([float(trck.PositionTooltipX1) * scale[0], float(trck.PositionTooltipY1) * scale[1],
                          float(trck.PositionTooltipZ1) * scale[2], float(trck.AngleX1), float(trck.AngleY1),
                          float(trck.AngleZ1)])
        reference = np.array([float(trck.PositionTooltipX2) * scale[0], float(trck.PositionTooltipY2) * scale[1],
                          float(trck.PositionTooltipZ2) * scale[2], float(trck.AngleX2), float(trck.AngleY2),
                          float(trck.AngleZ2)])

        if probe.all() and reference.all():
            coord = dynamic_reference(probe, reference)

    else:
        trck.Run()
        coord = np.array([float(trck.PositionTooltipX1) * scale[0], float(trck.PositionTooltipY1) * scale[1],
                          float(trck.PositionTooltipZ1) * scale[2], float(trck.AngleX1), float(trck.AngleY1),
                          float(trck.AngleZ1)])

    return coord


def PolhemusUSBCoord(trck, trck_id, ref_mode):
    endpoint = trck[0][(0, 0)][0]
    # Tried to write some settings to Polhemus in trackers.py while initializing the device.
    # TODO: Check if it's working properly.
    trck.write(0x02, "P")
    if trck_id == 2:
        scale = 10. * np.array([1., 1.0, -1.0])
    else:
        scale = 25.4 * np.array([1., 1.0, -1.0])
    coord = None

    if ref_mode:

        data = trck.read(endpoint.bEndpointAddress, 2 * endpoint.wMaxPacketSize)
        data = str2float(data.tostring())

        # six coordinates of first and second sensor: x, y, z and alfa, beta and gama
        # jump one element for reference to avoid the sensor ID returned by Polhemus
        probe = data[0] * scale[0], data[1] * scale[1], data[2] * scale[2], \
                data[3], data[4], data[5], data[6]
        reference = data[7] * scale[0], data[8] * scale[1], data[9] * scale[2], data[10], \
                    data[11], data[12], data[13]

        if probe.all() and reference.all():
            coord = dynamic_reference(probe, reference)

        return coord

    else:
        data = trck.read(endpoint.bEndpointAddress, endpoint.wMaxPacketSize)
        coord = str2float(data.tostring())

        coord = np.array((coord[0] * scale[0], coord[1] * scale[1], coord[2] * scale[2],
                          coord[3], coord[4], coord[5]))

        return coord


def PolhemusSerialCoord(trck_init, trck_id, ref_mode):
    # mudanca para fastrak - ref 1 tem somente x, y, z
    # aoflt -> 0:letter 1:x 2:y 3:z
    # this method is not optimized to work with all trackers, only with ISOTRAK
    # serial connection is obsolete, remove in future
    trck_init.write("P")
    lines = trck_init.readlines()

    coord = None

    if lines[0][0] != '0':
        print "The Polhemus is not connected!"
    else:
        for s in lines:
            if s[1] == '1':
                data = s
            elif s[1] == '2':
                data = s

        # single ref mode
        if not ref_mode:
            data = data.replace('-', ' -')
            data = [s for s in data.split()]
            j = 0
            while j == 0:
                try:
                    plh1 = [float(s) for s in data[1:len(data)]]
                    j = 1
                except:
                    print "error!!"

            coord = data[0:6]
    return coord


def DebugCoord(trk_init, trck_id, ref_mode):
    """
    Method to simulate a tracking device for debug and error check. Generate a random
    x, y, z, alfa, beta and gama coordinates in interval [1, 200[
    :param trk_init: tracker initialization instance
    :param ref_mode: flag for singular of dynamic reference
    :param trck_id: id of tracking device
    :return: six coordinates x, y, z, alfa, beta and gama
    """
    sleep(0.2)
    if ref_mode:
        probe = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200),
                          uniform(1, 200), uniform(1, 200), uniform(1, 200)])
        reference = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200),
                              uniform(1, 200), uniform(1, 200), uniform(1, 200)])

        coord = dynamic_reference(probe, reference)

    else:
        coord = np.array([uniform(1, 200), uniform(1, 200), uniform(1, 200),
                          uniform(1, 200), uniform(1, 200), uniform(1, 200)])

    return coord


def dynamic_reference(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])

    vet = probe[0:3] - reference[0:3]
    vet = np.mat(vet.reshape(3, 1))

    # Attitude Matrix given by Patriot Manual
    Mrot = np.mat([[cos(a) * cos(b), sin(b) * sin(g) * cos(a) - cos(g) * sin(a),
                       cos(a) * sin(b) * cos(g) + sin(a) * sin(g)],
                      [cos(b) * sin(a), sin(b) * sin(g) * sin(a) + cos(g) * cos(a),
                       cos(g) * sin(b) * sin(a) - sin(g) * cos(a)],
                      [-sin(b), sin(g) * cos(b), cos(b) * cos(g)]])

    coord_rot = Mrot.T * vet
    coord_rot = np.squeeze(np.asarray(coord_rot))

    return coord_rot[0], coord_rot[1], coord_rot[2], probe[3], probe[4], probe[5]


def str2float(data):
    """
    Converts string detected wth Polhemus device to float array of coordinates. THis method applies
    a correction for the minus sign in string that raises error while splitting the string into coordinates.
    :param data: string of coordinates read with Polhemus
    :return: six float coordinates x, y, z, alfa, beta and gama
    """

    count = 0
    for i, j in enumerate(data):
        if j == '-':
            data = data[:i + count] + ' ' + data[i + count:]
            count += 1

    data = [s for s in data.split()]
    data = [float(s) for s in data[1:len(data)]]

    return data