<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">"""

IN PROGRESS

Description:
This script was written to manually extract tissue thicknesses for the MULTIS in vivo and in vitro experimentation.
Anatomical and indentation trials are analyzed differently.
    Anatomical - All frames between the first and last pulse, arranges from minimum to maximum normal force
    Indentation - Frames from start of indentation to peak normal force of indentation

Getting started:
The user may either hard code the path to the subjectID folder or choose to open a file dialog to browse for the
folder.
The accepted trials are then loaded into a listbox.
Select the trial you would like to analyze by 'double-clicking' the name.
The first frame will appear for analysis.
Move the four red dots to the appropriate locations
    Superficial skin
    Skin/Fat boundary
    Fat/Muscle boundary
    Muscle/Bone boundary
To save results, press 'Enter' on the keyboard (Forces, moments, and thicknesses for that frame are saved to an xml
file)
The 'Next Image' button will take you to the next frame for analysis
When you are finished with analysis you can close the program one of two ways
    'Close Program' button in the lower right corner will close the windows and pretty print the newly created xml file
    The 'X' in the upper right corner will close all windows (without pretty printing the xml file)

    Original Author:
        Erica Morrill
        Department of Biomedical Engineering
        Lerner Research Institute
        Cleveland Clinic
        Cleveland, OH
        morrile2@ccf.org

"""

import matplotlib

matplotlib.use("TkAgg")
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import matplotlib.pyplot as plt
from matplotlib import patches
import os
import dicom
import SimpleITK as sitk
import tkFileDialog
from Tkinter import *
import Tkinter as tk
import tkMessageBox
import tdsmParserMultis
import numpy as np
import peakutils
import xml.etree.ElementTree as ET
import lxml.etree as etree


class DraggablePatch:
    lock = None  # only one can be animated at a time

    def __init__(self, obj):
        self.obj = obj
        self.press = None
        self.background = None

    def connect(self):
        """connect to all the events we need"""
        self.cidpress = self.obj.figure.canvas.mpl_connect('button_press_event', self.on_press)
        self.cidrelease = self.obj.figure.canvas.mpl_connect('button_release_event', self.on_release)
        self.cidmotion = self.obj.figure.canvas.mpl_connect('motion_notify_event', self.on_motion)

    def on_press(self, event):
        """on button press we will see if the mouse is over us and store some data"""
        if event.inaxes != self.obj.axes: return
        if DraggablePatch.lock is not None: return
        contains, attrd = self.obj.contains(event)
        if not contains: return

        x0, y0 = self.obj.center
        self.press = x0, y0, event.xdata, event.ydata

        DraggablePatch.lock = self

        # draw everything but the selected patch and store the pixel buffer
        canvas = self.obj.figure.canvas
        axes = self.obj.axes
        self.obj.set_animated(True)
        canvas.draw()
        self.background = canvas.copy_from_bbox(self.obj.axes.bbox)

        # now redraw just the patch
        axes.draw_artist(self.obj)

        # and blit just the redrawn area
        canvas.blit(axes.bbox)

    def on_motion(self, event):
        """on motion we will move the object if the mouse is over us"""
        if DraggablePatch.lock is not self:
            return
        if event.inaxes != self.obj.axes: return

        self.obj.center = (self.obj.center[0], event.ydata)

        canvas = self.obj.figure.canvas

        axes = self.obj.axes

        # restore the background region
        canvas.restore_region(self.background)

        # redraw just the current patch
        axes.draw_artist(self.obj)

        # blit just the redrawn area
        canvas.blit(axes.bbox)

    def on_release(self, event):
        """on release we reset the press data"""
        if DraggablePatch.lock is not self:
            return

        self.press = None
        DraggablePatch.lock = None

        # turn off the patch animation property and reset the background
        self.obj.set_animated(False)
        self.background = None

        # redraw the full figure
        self.obj.figure.canvas.draw()

    def disconnect(self):
        """disconnect all the stored connection ids"""
        self.obj.figure.canvas.mpl_disconnect(self.cidpress)
        self.obj.figure.canvas.mpl_disconnect(self.cidrelease)
        self.obj.figure.canvas.mpl_disconnect(self.cidmotion)

    def getlocation(self):
        return (self.obj.center[1])


def on_key(event, drs, mm_conv, TDMS_name, frame_num, data_zip):
    """Enter key was pressed, locations of each point (4) are found and sent to be calculated and saved to xml file"""
    if event.key == 'enter':
        coords = []
        for dr in drs:
            coord = dr.getlocation()
            coords.append(coord)
        S, F, M = calc_thicknesses(coords, mm_conv)  # Calculate the skin, fat and muscle thicknesses
        create_xml(S, F, M, TDMS_name, frame_num, data_zip)


def calc_thicknesses(coords, mm_conv):
    """Calculate thicknesses of skin, fat and muscle layers (in pixels) and convert to mm"""
    coords.sort()  # Sort the coordinates from lowest to highest to adjust for point switching

    skin = (coords[1] - coords[0]) * mm_conv
    fat = (coords[2] - coords[1]) * mm_conv
    muscle = (coords[3] - coords[2]) * mm_conv

    print(' ')
    print('*****Thickness Calculations*****')
    print('Skin_thickness = %f mm' % skin)
    print('Fat_thickness = %f mm' % fat)
    print('Muscle_thickness = %f mm' % muscle)
    print(' ')

    return skin, fat, muscle


def create_xml(S, F, M, TDMS_name, f_num, data_zip):
    """Create xml file with forces, moments, and thicknesses"""
    global mm
    # print('Saved data in XML')

    #Extract data for specific frame (i.e. 'mm')
    data = data_zip[mm]
    fx = data[1]
    fy = data[2]
    fz = data[3]
    mx = data[4]
    my = data[5]
    mz = data[6]

    #Define path for the analysis folder
    split_name = os.path.split(TDMS_name)
    analysis_path = os.path.dirname(os.path.dirname(TDMS_name)) + '/Analysis'

    # Check for Analysis directory
    if not os.path.isdir(analysis_path):
        os.makedirs(analysis_path)

    # Name of the xml file
    xml_name = os.path.dirname(os.path.dirname(TDMS_name)) + '/Analysis/' + split_name[1][4:25] + '.xml'

    # Check if an xml file exists for this image. If TRUE, edit the values and if FALSE, create a new one.
    if os.path.exists(xml_name):
        # print('Edit File')
        doc = ET.parse(xml_name)
        root = doc.getroot()
        loc_lst = []
        for loc in root.findall('Location'):
            loc_lst.append(loc.find("Name").text)
        tail = split_name[1][16:25] + '-Frame_' + str(f_num)

        if tail in loc_lst:
            # Over-write thicknesses: Forces and moments should be the same since the data for this frame has already
            # been recorded
            loc = root.find('Location[Name="%s"]' % tail)
            Thick = root.find("Thickness")
            Thick.find("Skin").text = str(S)
            Thick.find("Fat").text = str(F)
            Thick.find("Muscle").text = str(M)

            tree = ET.ElementTree(root)
            tree.write(xml_name, xml_declaration=True)

        else:
            #Add new child, because information for this frame is not yet in the xml file
            loc = ET.SubElement(root, 'Location')
            ET.SubElement(loc, "Name").text = tail

            Forces = ET.SubElement(loc, "Forces")
            ET.SubElement(Forces, "Fx", units="N").text = str(fx)
            ET.SubElement(Forces, "Fy", units="N").text = str(fy)
            ET.SubElement(Forces, "Fz", units="N").text = str(fz)

            Moments = ET.SubElement(loc, "Moments")
            ET.SubElement(Moments, "Mx", units="Nm").text = str(mx)
            ET.SubElement(Moments, "My", units="Nm").text = str(my)
            ET.SubElement(Moments, "Mz", units="Nm").text = str(mz)

            Thick = ET.SubElement(loc, "Thickness")
            ET.SubElement(Thick, "Skin", units="mm").text = str(S)
            ET.SubElement(Thick, "Fat", units='mm').text = str(F)
            ET.SubElement(Thick, "Muscle", units='mm').text = str(M)

            tree = ET.ElementTree(root)
            tree.write(xml_name, xml_declaration=True)

    else:
        # print('Make New File')

        tail = split_name[1][16:25] + '-Frame_' + str(f_num)

        subID = os.path.split(os.path.dirname(os.path.dirname(TDMS_name)))
        root = ET.Element(subID[1])
        loc = ET.SubElement(root, 'Location')
        ET.SubElement(loc, "Name").text = tail

        Forces = ET.SubElement(loc, "Forces")
        ET.SubElement(Forces, "Fx", units="N").text = str(fx)
        ET.SubElement(Forces, "Fy", units="N").text = str(fy)
        ET.SubElement(Forces, "Fz", units="N").text = str(fz)

        Moments = ET.SubElement(loc, "Moments")
        ET.SubElement(Moments, "Mx", units="Nm").text = str(mx)
        ET.SubElement(Moments, "My", units="Nm").text = str(my)
        ET.SubElement(Moments, "Mz", units="Nm").text = str(mz)

        Thick = ET.SubElement(root, "Thickness")
        ET.SubElement(Thick, "Skin", units="mm").text = str(S)
        ET.SubElement(Thick, "Fat", units='mm').text = str(F)
        ET.SubElement(Thick, "Muscle", units='mm').text = str(M)

        tree = ET.ElementTree(root)
        tree.write(xml_name, xml_declaration=True)


def make3digits(input):
    """Convert all trial numbers to three digits (if they are not already)"""
    if len(input) == 1:
        output = '00' + str(input)
    elif len(input) == 2:
        output = '0' + str(input)
    else:
        output = str(input)
    return output


def createAverageFit(Fx, avgThres):
    """Filter the tdms normal force data"""
    avgeragelist = []
    for items in range(len(Fx)):
        if Fx[items] != Fx[items - avgThres]:
            num2avg = Fx[items:(items + avgThres)]
            avgeragelist.append(np.average(num2avg))
        else:
            continue

    avgeragelist = [avgeragelist[0]] * (avgThres / 2) + avgeragelist[0:-(avgThres) / 2]
    return avgeragelist


def findFrame(initial_time, dT, frameTimeVector):
    """Find the frame corresponding to the specified tdms time and return the adjusted tdms time to match the
    selected frame"""
    adjusted_time = dT + initial_time
    for frames in range(len(frameTimeVector)):
        frames += 1
        frame_time = sum(frameTimeVector[0:frames])
        if adjusted_time &lt;= frame_time:
            timeDiff_up = frame_time - adjusted_time
            timeDiff_low = adjusted_time - sum(frameTimeVector[0:frames - 1])
            if timeDiff_up &lt; timeDiff_low:
                frame_frame = frames
                readjusted_time_tdms = frame_time - dT
            else:
                frame_frame = frames - 1
                readjusted_time_tdms = sum(frameTimeVector[0:frames - 1]) - dT
            break

    return frame_frame, int(readjusted_time_tdms)


def getFiles():
    """Get accepted trials for the selected subject, this includes the Ultrasound dicom images, Data tdms files,
    and the TimeSynchronization txt file. The program can either open a File Dialog box to search for the directory,
    or the user can hard code the directory"""
    # #Use File Dialog to open Subject Directory
    # root = Tk()
    # directoryname = tkFileDialog.askdirectory(title = 'Please select a SubjectID folder to be analyzed')
    # root.destroy()

    # Hard code Subject Directory
    directoryname = '/home/morrile2/Documents/MULTIS_test/MULTIS001-1'

    timeFile = open(directoryname + '/TimeSynchronization.txt', 'r')
    lines = timeFile.readlines()

    trial = []
    delta_t = []

    for i in lines:
        l = i.split(' ')
        trial.append(l[0])
        delta_t.append(l[1])
    timeFile.close()

    for ii in range(len(trial)):
        trial[ii] = make3digits(trial[ii])

    time_adj = zip(trial, delta_t)

    lstFilesDCM = []  # create an empty list for DICOM files
    lstFilesTDMS = []  # create an empty list for TDMS files

    for dirName, subdirList, fileList in os.walk(directoryname):
        for filename in fileList:
            if filename[0:3] in trial:
                if ".ima" in filename.lower():  # check whether the file's DICOM
                    lstFilesDCM.append(os.path.join(dirName, filename))
                elif ".tdms" in filename.lower():  # check whether the file's TDMS
                    lstFilesTDMS.append(os.path.join(dirName, filename))

    # Sort variables by trial number
    lstFilesTDMS.sort()
    lstFilesDCM.sort()
    time_adj.sort()

    return lstFilesDCM, lstFilesTDMS, time_adj


def getFrames(file_DCM, file_TDMS, time_adj):
    """Get the frames to be analyzed and save the data for those frames. Different for indentation and anatomical
    trials. Indentation contains frames that start at indentation and go through the peak force of the indentation,
    while anatomical analyzes all frames between start and end pulses from minimum to maximum force"""

    # Extract force information from TDMS file
    data = tdsmParserMultis.parseTDMSfile(file_TDMS)

    Fx = np.array(data[u'State.6-DOF Load'][u'6-DOF Load Fx'])
    Fy = np.array(data[u'State.6-DOF Load'][u'6-DOF Load Fy'])
    Fz = np.array(data[u'State.6-DOF Load'][u'6-DOF Load Fz'])
    Mx = np.array(data[u'State.6-DOF Load'][u'6-DOF Load Mx'])
    My = np.array(data[u'State.6-DOF Load'][u'6-DOF Load My'])
    Mz = np.array(data[u'State.6-DOF Load'][u'6-DOF Load Mz'])

    pulse = np.array(data[u'Sensor.Run Number Pulse Train'][u'Run Number Pulse Train'])
    pulse = pulse[:]

    Peaks = peakutils.indexes(pulse, thres=0.5 * max(pulse), min_dist=100)

    dT = float(time_adj[1])

    # Read metadata from the dicom image
    RefDs = dicom.read_file(file_DCM)
    frameTimeVector = RefDs.FrameTimeVector
    seq = RefDs.SequenceofUltrasoundRegions
    conv_fact = seq[0].PhysicalDeltaY * 10  # Extract conversion factor from dicom metadata (Original is in cm, convert
    # to mm)

    if file_TDMS[-8:-7] == 'I':
        indentation = True
        anatomical = False
    elif file_TDMS[-8:-7] == 'A':
        indentation = False
        anatomical = True

    force_list = list(Fx)

    time = np.arange(len(Fx))
    data_i = zip(time, Fx, Fy, Fz, Mx, My, Mz)

    if indentation:
        # Indentation frames from start of indentation to peak force during indentation
        # Find indentation start time in TDMS timeline , denoted by tdms b/c used for tdms
        avg = createAverageFit(Fx, 300)
        index_range = 200
        for items in range(len(avg[0:-index_range])):
            if abs(avg[items + index_range] - avg[items]) &gt; 1.5:
                time_start_tdms = items + index_range / 2
                break

        time_max_tdms = force_list.index(min(force_list))  # note max and min reversed as Fx is negative force

        start_frame, start_frame_time_tdms = findFrame(time_start_tdms, dT, frameTimeVector)
        max_frame, max_frame_time_tdms = findFrame(time_max_tdms, dT, frameTimeVector)

        frame_lst_i = [start_frame]
        data_i_final = [data_i[start_frame_time_tdms]]

        for t in data_i[start_frame_time_tdms + 1:max_frame_time_tdms]:
            inc_frame, inc_frame_time_tdms = findFrame(t[0], dT, frameTimeVector)
            if inc_frame not in frame_lst_i:
                frame_lst_i.append(inc_frame)
                data_i_final.append(data_i[inc_frame_time_tdms])

        frames = frame_lst_i
        DATA = data_i_final

    elif anatomical:
        # Anatomical frames from minimum to maximum normal force (Fx)
        time_preIndent_tdms = Peaks[0]  # 230 ms is location first pulse.
        time_postIndent_tdms = Peaks[-1]

        # Sort the data by the normal force (Fx) to get anatomical frame list from minimum to maximum
        data_a = zip(time, Fx, Fy, Fz, Mx, My, Mz)
        data_a.sort(key=lambda t: abs(t[1]))

        frame_lst_a = []
        data_a_sort = []

        for t in data_a:
            min_frame, min_frame_time_tdms = findFrame(t[0], dT, frameTimeVector)
            if min_frame not in frame_lst_a:
                if min_frame_time_tdms &gt; time_preIndent_tdms and min_frame_time_tdms &lt; time_postIndent_tdms:
                    frame_lst_a.append(min_frame)
                    data_a_sort.append(data_i[min_frame_time_tdms])

        final_a_sort = zip(frame_lst_a, data_a_sort)
        final_a_sort.sort(key=lambda t: abs(t[1][1]))
        frames = []
        DATA = []
        for d in final_a_sort:
            frames.append(d[0])
            DATA.append(d[1])

    return frames, conv_fact, DATA


def next_image(button1):
    """Gets next image for analysis - Triggered by 'Next Image' button"""
    global location, mm, selectedFrames, root_plot
    mm += 1
    if mm &lt; len(selectedFrames):
        root_plot.destroy()
        execute(location)
    elif mm == len(selectedFrames) - 1:
        button1.config(state="disabled")


def end_program(TDMS_name):
    """Checks to see if there was an xml file created. If so, it is pretty printed to the filename and the program
    closes all windows. Triggered by the 'Close Program' button"""
    global app, root_plot

    split_name = os.path.split(TDMS_name)
    xml_name = os.path.dirname(os.path.dirname(TDMS_name)) + '/Analysis/' + split_name[1][4:25] + '.xml'

    if os.path.exists(xml_name):
        prettyPrintXml(xml_name)

    root_plot.destroy()
    app.quit()


def prettyPrintXml(xmlFilePathToPrettyPrint):
    """Pretty print the xml file after all frames have been analyzed"""
    assert xmlFilePathToPrettyPrint is not None
    parser = etree.XMLParser(resolve_entities=False, strip_cdata=False)
    document = etree.parse(xmlFilePathToPrettyPrint, parser)
    document.write(xmlFilePathToPrettyPrint, pretty_print=True, encoding='utf-8')


def on_exit():
    """Closes all windows if the user selects okay, does not pretty print the xml file"""
    global app, root_plot
    if tkMessageBox.askokcancel("Quit", "Do you want to quit without pretty-printing xml file?"):
        root_plot.destroy()
        app.quit()


def showPlot(DCM_name, img_idx, conv_fact, TDMS_name, data_zip):
    """Read and show the appropriate frame of the dicom image in tkinter GUI"""
    global img_all, mm, root_plot

    # Read dicom image
    if mm == 0:
        reader = sitk.ImageFileReader()
        reader.SetFileName(DCM_name)
        img_all = reader.Execute()

        # Plot the figure
        figsize = (18, 12)
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(111)
    else:
        fig = plt.gcf()
        fig.clear()
        ax = fig.add_subplot(111)

    # Read in the appropriate dicom frame and set up plot title
    img = img_all[:, :, img_idx]
    nda = sitk.GetArrayFromImage(img)
    plt.imshow(nda[:, :, 0], cmap="gray")

    split_name = os.path.split(TDMS_name)
    plt.title(split_name[1][0:25] + '  :  Frame = ' + str(img_idx) + '\n\n' + str(mm + 1) + '  of  ' + str(
        len(data_zip)) + '\n')

    # Create four circles used to indicate tissue interfaces (Top of skin, skin/fat, fat/muscle, bottom of muscle)
    r = 7
    circs = [matplotlib.patches.Circle((512, 60), radius=r, alpha=0.5, fc='r'),
             matplotlib.patches.Circle((512, 90), radius=r, alpha=0.5, fc='r'),
             matplotlib.patches.Circle((512, 160), radius=r, alpha=0.5, facecolor='r'),
             matplotlib.patches.Circle((512, 360), radius=r, alpha=0.5, facecolor='r')]
    drs = []

    # Create main GUI
    root_plot = tk.Tk()
    root_plot.title("Drag points to tissue boundaries and hit 'Enter' to save")
    canvas_plot = FigureCanvasTkAgg(fig, master=root_plot)
    canvas_plot.draw()
    canvas_plot.get_tk_widget().grid(row=0, column=0)

    # Next image, grabs the next frame in the sequence to analyze for this dicom image
    button1 = tk.Button(root_plot, text="Next Image", command=lambda: next_image(button1)).grid(row=2, column=0)
    # Close Program, ends the script, closes all windows and writes the xml file using pretty print
    tk.Button(root_plot, text="Close Program", command=lambda: end_program(TDMS_name)).grid(row=2, column=0, sticky="E")
    # When the main window is closed, the user is asked if they want to close without executing pretty print for the
    # xml file
    root_plot.protocol("WM_DELETE_WINDOW", on_exit)

    # Plot the four circles and make them draggable using the DraggablePatch Class
    for circ in circs:
        ax.add_patch(circ)
        dr = DraggablePatch(circ)
        dr.connect()
        drs.append(dr)

    # Press enter when you are happy with circle locations to calculate and record skin, fat, and muscle thicknesses
    cid = fig.canvas.mpl_connect('key_press_event',
                                 lambda event: on_key(event, drs, conv_fact, TDMS_name, img_idx, data_zip))


def execute(kk):
    """Initiate the main application to show the Ultrasound images and interactive GUI for the thickness measurements"""
    global lstFilesDCM, lstFilesTDMS, time_adj, mm, location, selectedFrames, conv_fact, data, app

    location = kk

    if mm == 0:
        app.destroy()  # Closes list window because the file selection is complete
        selectedFrames, conv_fact, data = getFrames(lstFilesDCM[location], lstFilesTDMS[location], time_adj[location])

    selectedFrames[0]
    if mm &lt; len(selectedFrames):
        showPlot(lstFilesDCM[location], selectedFrames[mm], conv_fact, lstFilesTDMS[location], data)


class FileSelectionApp(tk.Tk):
    """Application to display all of the accepted trials in a list"""

    def __init__(self, *args, **kwargs):
        global lstFilesTDMS
        tk.Tk.__init__(self, *args, **kwargs)
        self.title('Select Location for Analysis')
        lb = tk.Listbox(self)
        sb = tk.Scrollbar(self, orient=VERTICAL, command=lb.yview)
        lb.config(yscrollcommand=sb.set)

        lb.grid(column=0, row=0, sticky=(N, W, E, S))
        sb.grid(column=1, row=0, sticky=(N, S))
        lb.config(width=40, height=30)

        for item in lstFilesTDMS:
            lb.insert("end", item[-30:-5])

        lb.bind("&lt;Double-Button-1&gt;", self.OnDouble)

    def OnDouble(self, event):
        """When a location name is specified using a double-click, the execute function begins analysis"""
        widget = event.widget
        selection = widget.curselection()
        value = map(int, selection)
        execute(value[0])

    def yview(self, *args):
        apply(self.yview, args)


if __name__ == "__main__":
    # Grabs subject folder and initiates the beginning of the application (i.e. List of accepted trials)
    global lstFilesTDMS, lstFilesDCM, time_adj, mm, location, app
    mm = 0
    lstFilesDCM, lstFilesTDMS, time_adj = getFiles()  # Get all the associated accepted files for subject
    app = FileSelectionApp()
    app.mainloop()
</pre></body></html>