import numpy as np import sys from tvtk.api import tvtk from tvtk.api import colors from lxml import etree as et import os import subprocess import warnings warnings.simplefilter(action='ignore', category=FutureWarning) sys.path.append('/home/landisb/PycharmProjects/PythonScripts/Project') import myVTK import RunPygem def CreateNestedDictionary(): """Create a nested dictionary structure for the ultrasound data this iterates "nicely" for these purposes as for place, slice in nestedDict.items(): for orient, strata in slice.items(): for layer, point in strata.items():""" from copy import deepcopy import numpy as np orient = { 'bone' : np.zeros(3), # 'muscle': np.zeros(3), # 'fat' : np.zeros(3), 'skin' : np.zeros(3) } place = { 'Anterior' : deepcopy(orient), 'Posterior': deepcopy(orient), 'Medial' : deepcopy(orient), 'Lateral' : deepcopy(orient) } nestedDict={'Proximal' : deepcopy(place), 'Central' : deepcopy(place), 'Distal' : deepcopy(place) } return nestedDict def GetSubjectData(filename): """Reads in-vivo subject file MULTIS00X-1.xml returns subjectData in a dictionary""" subjectData = {} tree = et.ElementTree(file=filename) root = tree.getroot() subject = root.find('Subject_Data') measurements = subject.find('Anatomical_Measurements') limbs = [] limbs.append(measurements.find('Lower_Leg') ) limbs.append(measurements.find('Lower_Arm')) limbs.append(measurements.find('Upper_Leg')) limbs.append(measurements.find('Upper_Arm')) # for limb in limbs: subjectData[limb.tag]={} for cluster in limb: children = list(cluster) descriptor, value = children[0].text, children[1].text descriptor = descriptor.replace(' ','') annoyingDescriptions =('Hip', 'Knee', 'Elbow', 'Shoulder') for annoy in annoyingDescriptions: descriptor=descriptor.replace(annoy,'') subjectData[limb.tag][descriptor] = float(value) return subjectData def GetCadaverData(filename): """Reads in-vitro cadaver subject file CMULTIS00X-1.xml this file is structured slightly differently than in-vivo returns subjectData in a dictionary""" subjectData = {} tree = et.ElementTree(file=filename) root = tree.getroot() subject = root.find('Subject_Data') limbs = ['Arm', 'Leg'] segments = ['Upper', 'Lower'] measurements = [] for limb in limbs: temp = subject.find('{}_-_Anatomical_Measurements'.format(limb)) for segment in segments: measurements.append(temp.find('{}_{}'.format(segment,limb))) for measure in measurements: subjectData[measure.tag]={} for cluster in measure: children = list(cluster) descriptor, value = children[0].text, children[1].text descriptor = descriptor.replace(' ','') annoyingDescriptions =('Hip', 'Knee', 'Elbow', 'Shoulder') for annoy in annoyingDescriptions: descriptor=descriptor.replace(annoy,'') subjectData[measure.tag][descriptor] = float(value) return subjectData def ReadUltrasoundThickness(limbSegment, xmlFile): """Read the xml to find thicknesses""" import numpy as np import os import re from lxml import etree as et directory = os.path.dirname(xmlFile) doc = et.parse(xmlFile) root = doc.getroot() bodyPart = root.find(limbSegment.replace('_', '')) AllThicknesses = CreateNestedDictionary() for loc in list(bodyPart): place, orient = re.findall('[A-Z][a-z]*', loc.tag) # Split on capitalization thicknessFile = loc.get('Anatomical') subDoc = et.parse(os.path.join(directory, thicknessFile)) subRoot = subDoc.getroot() subject = subRoot.find("Subject").find("Source") Thickness = {} for frame in list(subject): readings = frame.find('Thickness') for value in list(readings): try: Thickness[value.tag.lower()].append(float(value.text)) except KeyError: Thickness[value.tag.lower()] = [float(value.text)] for t, v in Thickness.iteritems(): if np.all(np.isnan(v)): # no data, use CT thickness evaluation instead # except I don't have this for in-vivo data, only my test case... pass # v = [] # empty the list as it is only NaNs # # note the follwing code is similar but NOT the same as the subDoc section above. # subDoc = et.parse(os.path.join(directory, '../CTManual_zProbe',thicknessFile)) # adjust path to CT thickness directory # subRoot = subDoc.getroot() # subject = subRoot.find("Subject").find("Source") # for frame in list(subject): # value = frame.find('Thickness').find(t.capitalize()) # v.append(float(value.text)) elif np.any(np.isnan(v)): v = v[~np.isnan(v)] # just use the good values Thickness[t] = sum(v) / len(v) AllThicknesses[place][orient] = Thickness return AllThicknesses def ReadInVitroPositions(limbSegment, xmlFile): """Read the xml to find the in vitro ultrasound positions""" import re from lxml import etree as et basePoints = CreateNestedDictionary() baseDirections = CreateNestedDictionary() newFile = xmlFile.replace('/','\\') print newFile, os.path.isfile(newFile) doc = et.parse(xmlFile) root = doc.getroot() bodyPart = root.find(limbSegment.replace('_','')) for loc in list(bodyPart): place = loc.tag place, orient = re.findall('[A-Z][a-z]*', place) # Split on capitalization pos = loc.find("Anatomical").find("USPosition") x = float( pos.find("x").get('value') ) * 1000 y = float( pos.find("y").get('value') ) * 1000 z = float( pos.find("z").get('value') ) * 1000 roll = float( pos.find("roll" ).get('value') ) pitch= float( pos.find("pitch").get('value') ) yaw = float( pos.find("yaw" ).get('value') ) basePoints[place][orient]['skin'] = [x, y, z] baseDirections[place][orient] = findDirectionFromRPY(roll, pitch, yaw) return basePoints, baseDirections def findDirectionFromRPY(roll, pitch, yaw, reverse=False): from numpy import sin, cos, dot, array r = roll p = pitch y = yaw transformMatrix = array( [[ cos(p)*cos(y), sin(r)*sin(p)*cos(y) - cos(r)*sin(y), cos(r)*sin(p)*cos(y) + sin(r)*sin(y)], [ cos(p)*sin(y), sin(r)*sin(p)*sin(y) + cos(r)*cos(y), cos(r)*sin(p)*sin(y) - sin(r)*cos(y)], [-sin(p), sin(r)*cos(p), cos(r)*cos(p) ] ]) if reverse: dirVector = dot(transformMatrix, [0, 0,-1]) else: dirVector = dot(transformMatrix, [0, 0, 1]) # ultrasound probe points in +Z dir return dirVector def SimpleMorphing(renderer, subjectDir, cadaverDir, meshes, limb, activeMorphing=True, zPlane=False): """Morph a model (cadaver) into a shape matching an in-vivo subject""" # draw meshes myVTK.MakeStl(ren, meshes['bone'], color=colors.blanched_almond, opacity=0.1) myVTK.MakeStl(ren, meshes['skin'], color=colors.cyan, opacity=0.1) # Create handle for morphed meshes morphs = {layer: mesh.replace('.','_morphed_{}.'.format(subjectDir.rsplit(os.sep,1)[1])) for layer,mesh in meshes.items() if layer is not 'bone'} # # Find files from positions relative to given directories subjectFile = os.path.join(subjectDir, 'Configuration', '{}.xml'.format(subjectDir.rsplit(os.sep, 1)[1])) subjectUltrasound = os.path.join(subjectDir, 'TissueThickness', 'UltrasoundManual', '{}_TA_inclusion.xml'.format(subjectDir.rsplit(os.sep, 1)[1])) cadaverFile = os.path.join(cadaverDir, 'Configuration', '{}.xml'.format(cadaverDir.rsplit(os.sep, 1)[1])) cadaverUltrasound = os.path.join(cadaverDir, 'TissueThickness', 'UltrasoundManual', '{}_TA_inclusion.xml'.format(cadaverDir.rsplit(os.sep, 1)[1])) inVitroPositionsFile = os.path.join(cadaverDir, 'Registration', '{}_{}_US_CT.xml'.format(cadaverDir.rsplit(os.sep, 1)[1], ''.join(L for L in limb if L.isupper()))) # # Read cadaver subject file to get in-vitro experimentally meased values cadaverMeasurements = GetCadaverData(cadaverFile) cadaverLength = cadaverMeasurements[limb]['Length'] * 10 cadaverCircum = { 'Central' : cadaverMeasurements[limb][ 'CentralCircumference']*10, 'Distal' : cadaverMeasurements[limb][ 'DistalCircumference']*10, 'Proximal': cadaverMeasurements[limb]['ProximalCircumference']*10 } # Read subject file xml to find in-vivo experimentally measured sizes subjectMeasurements = GetSubjectData(subjectFile) subjectLength = subjectMeasurements[limb]['Length']*10 subjectCircum = { 'Central' : subjectMeasurements[limb][ 'CentralCircumference']*10, 'Distal' : subjectMeasurements[limb][ 'DistalCircumference']*10, 'Proximal': subjectMeasurements[limb]['ProximalCircumference']*10 } # # distance from central plane, file is distance from hip...so math is done now subjectDistances= { 'Hip' : subjectMeasurements[limb]['landmarktoCentralCircumference']*10, 'Proximal': subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectMeasurements[limb]['landmarktoProximalCircumference']*10, 'Central' : 0.0, 'Distal' : subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectMeasurements[limb]['landmarktoDistalCircumference']*10, 'Knee' : subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectLength } # # Read ultrasound thicknesses produces layer thickness cadaverThicknesses = ReadUltrasoundThickness(limb, cadaverUltrasound) subjectThickness = ReadUltrasoundThickness(limb, subjectUltrasound) # sum layers for total thickness cadaverThicknesses= {place:{orient:sum(strata.values()) for orient, strata in slice.items() } for place,slice in cadaverThicknesses.iteritems()} subjectThickness = {place:{orient:sum(strata.values()) for orient, strata in slice.items() } for place,slice in subjectThickness.iteritems()} # # Read probe positions probePoints, probeDirections = ReadInVitroPositions(limb, inVitroPositionsFile) probeDirections = {place:{orient:strata for orient, strata in slice.items()} for place, slice in probeDirections.iteritems()} probePoints = {place:{orient:[point for layer,point in strata.items()if layer is 'skin'][0] for orient,strata in slice.items()}for place,slice in probePoints.items()} # if zPlane: # if true forces probePoints into the same z plane for place, slice in probePoints.iteritems(): for orient, strata in slice.items(): strata[2] = probePoints[place]['Anterior'][2] # # Find control point for morph controlPoints = CreateNestedDictionary() modelThicknesses = CreateNestedDictionary() norms = { 'Anterior' : np.array([ 0,-1, 0]), 'Lateral' : np.array([-1, 0, 0]), 'Medial' : np.array([ 1, 0, 0]), 'Posterior' : np.array([ 0, 1, 0]) } for place, slice in probePoints.iteritems(): # if place is 'Central': for orient, strata in slice.items(): p1 = probePoints[place][orient] p2 = p1 + norms[orient] * cadaverThicknesses[place][orient] p2 = p1 + probeDirections[place][orient] * cadaverThicknesses[place][orient] # controlPoints[place][orient]['bone'] = myVTK.FindIntersectionOfStlAndLine(ren, meshes['bone'], p1, p2, colors.red , visualize=True, embellish=True, extend=True) controlPoints[place][orient]['skin'] = myVTK.FindIntersectionOfStlAndLine(ren, meshes['skin'], p1, p2, colors.green, visualize=False, embellish=True, extend=True) modelThicknesses[place][orient] = myVTK.dist(controlPoints[place][orient]['bone'], controlPoints[place][orient]['skin']) # # myVTK.Show(ren, mouseHighlightActor=False) # morphedPoints = CreateNestedDictionary() for place, slice in controlPoints.items(): for orient, strata in slice.items(): for layer, point in strata.items(): if layer is 'bone': morphedPoints[place][orient][layer] = point else: if np.isnan(subjectThickness[place][orient]): morphedPoints[place][orient][layer] = controlPoints[place][orient][layer] else: morphedPoints[place][orient][layer] = controlPoints[place][orient]['bone']-norms[orient]*subjectThickness[place][orient] print place, orient, layer, ' ',morphedPoints[place][orient][layer] # controlArray = np.array( [point for place,slice in controlPoints.items() for orient,strata in slice.items() for layer,point in strata.items() if layer is not 'bone'] ) morphedArray = np.array( [point for place,slice in morphedPoints.items() for orient,strata in slice.items() for layer,point in strata.items() if layer is not 'bone'] ) # # Morph meshes with subprocess to speed things up. Morphing of seperate layers will run concurrently. if activeMorphing: # This was used for testing purposes, it is left in because the user may want the output again, without recalculating the morph print 'Morphing...', processes = [] for layer, mesh in meshes.items(): if layer is not 'bone': RunPygem.WritePygemParameterFile('{}.dat'.format(layer), controlArray,morphedArray) processes.append( subprocess.Popen(['python', 'RunPygem.py', mesh, '{}.dat'.format(layer), morphs[layer]]) ) print 'this takes a while...', for proc in processes: proc.wait() print 'finished.' # myVTK.MakeStl(ren, morphs['skin'], color=colors.orange, opacity=0.1) # morphThicknesses = CreateNestedDictionary() for place, slice in morphedPoints.iteritems(): if place is 'Central': for orient, strata in slice.items(): p1 = morphedPoints[place][orient]['skin'] p2 = p1 + norms[orient]*(modelThicknesses[place][orient]+subjectThickness[place][orient])/2 # intersectBone = myVTK.FindIntersectionOfStlAndLine(ren, meshes['bone'], p1, p2, colors.red , embellish=True, extend=True) # no longer morphing bone intersectSkin = myVTK.FindIntersectionOfStlAndLine(ren, morphs['skin'], p1, p2, colors.blue, embellish=True, extend=True, size=1) morphThicknesses[place][orient] = myVTK.dist(intersectBone, intersectSkin) diff = morphThicknesses[place][orient]-subjectThickness[place][orient] # myVTK.Show(ren, mouseHighlightActor=False) # planeNormal = (0,0,1) planeSize = 250 inVitroColors = {'Proximal':colors.pink, 'Central':colors.green_pale, 'Distal':colors.cyan, 'Hip':colors.slate_grey, 'Knee':colors.yellow} inVivo_Colors = {'Proximal':colors.red, 'Central':colors.green_dark, 'Distal':colors.blue, 'Hip':colors.black, 'Knee':colors.orange} # modelCircum = {} morphCircum = {} for place, slice in controlPoints.items(): if place is 'Central': cP = controlPoints[place]['Anterior']['skin'] # makes following lines somewhat more readable mP = morphedPoints[place]['Anterior']['skin'] # myVTK.DrawPlane(ren, cP, planeNormal, size=planeSize, opacity=0.05, color=inVitroColors[place]) # for visualization only plane = tvtk.Plane(origin=cP, normal=planeNormal) # Create the plane as a function, that will be used to cut the stl. modelCircum[place],c = myVTK.MakeStlCut(ren, meshes['skin'], plane, inVitroColors[place], opacity=0) # # myVTK.DrawPlane(ren, mP, normal, size=planeSize, opacity=0.05, color=inVivo_Colors[place]) # for visualization only plane = tvtk.Plane(origin=mP, normal=planeNormal) # Create the plane as a function, that will be used to cut the stl. morphCircum[place],c = myVTK.MakeStlCut(ren, morphs['skin'], plane, inVivo_Colors[place], opacity=0) # MBoneCircum[place],c = myVTK.MakeStlCut(ren, morphs['bone'], plane, colors.black, opacity=1) # print '{:8s} Bone Morph circum {:6.2f}, Bone circum {:6.2f}'.format(place, MBoneCircum[place], bone_Circum[place]) print place # print ' Model circum {: 6.2f} thickness {: 6.2f}'.format(modelCircum[place], modelThicknesses[place]) print ' Morph circum {: 6.2f}'.format(morphCircum[place]), print ' Thickness:', print ' Anterior {: 6.2f}'.format(morphThicknesses[place]['Anterior']), print ' Posterior {: 6.2f}'.format(morphThicknesses[place]['Posterior']), print ' Lateral {: 6.2f}'.format(morphThicknesses[place]['Lateral']), print ' Medial {: 6.2f}'.format(morphThicknesses[place]['Medial']) # print ' Subj circum {: 6.2f}'.format(subjectCircum[place]), print ' Thickness:', print ' Anterior {: 6.2f}'.format(subjectThickness[place]['Anterior'] ), print ' Posterior {: 6.2f}'.format(subjectThickness[place]['Posterior']), print ' Lateral {: 6.2f}'.format(subjectThickness[place]['Lateral']), print ' Medial {: 6.2f}'.format(subjectThickness[place]['Medial']) # print ' Cadaver circum {: 6.2f}'.format(cadaverCircum[place]), print ' Thickness:', print ' Anterior {: 6.2f}'.format(cadaverThicknesses[place]['Anterior'] ), print ' Posterior {: 6.2f}'.format(cadaverThicknesses[place]['Posterior']), print ' Lateral {: 6.2f}'.format(cadaverThicknesses[place]['Lateral']), print ' Medial {: 6.2f}'.format(cadaverThicknesses[place]['Medial']) # print ' Model circum {: 6.2f}'.format(modelCircum[place]), print ' Thickness:', print ' Anterior {: 6.2f}'.format(modelThicknesses[place]['Anterior'] ), print ' Posterior {: 6.2f}'.format(modelThicknesses[place]['Posterior']), print ' Lateral {: 6.2f}'.format(modelThicknesses[place]['Lateral']), print ' Medial {: 6.2f}'.format(modelThicknesses[place]['Medial']) if __name__ == '__main__': # gender = 'female' limb = 'Lower_Leg' # choose from Upper_Arm, Upper_Leg, Lower_Arm, Lower_Leg # # This is a coded reminder of the model and invivo parings if gender is 'female': subjectDir = 'C:\Users\landisb\Workspace\InVivoSubjects\MULTIS037-1' cadaverDir = 'C:\Users\landisb\Workspace\CadaverModels\CMULTIS008-1' elif gender is 'male': subjectDir = 'C:\Users\landisb\Workspace\InVivoSubjects\MULTIS094-1' cadaverDir = 'C:\Users\landisb\Workspace\CadaverModels\CMULTIS006-1' # # Adjust meshes to match cadaver and limb segment choice meshes = { 'bone': r'C:\Users\landisb\Downloads\Operation Multis-STL_Files\STL_Files\Lumped\008LL\meshed_bone_quad.stl', 'skin': r'C:\Users\landisb\Downloads\Operation Multis-STL_Files\STL_Files\Lumped\008LL\008LL_meshed_lump.stl' } # ren = tvtk.Renderer(background=colors.white) SimpleMorphing(ren, subjectDir, cadaverDir, meshes, limb, zPlane=False, activeMorphing=True) print 'complete' myVTK.Show(ren, mouseHighlightActor=False)