1 import numpy as np
   2 import sys
   3 from tvtk.api import tvtk
   4 from tvtk.api import colors
   5 from lxml import etree as et
   6 import os
   7 import subprocess
   8 import warnings
   9 warnings.simplefilter(action='ignore', category=FutureWarning)
  10 
  11 sys.path.append('/home/landisb/PycharmProjects/PythonScripts/Project')
  12 import myVTK
  13 import RunPygem
  14 
  15 
  16 def CreateNestedDictionary():
  17 	"""Create a nested dictionary structure for the ultrasound data
  18 		this iterates "nicely" for these purposes as
  19 		for place, slice in nestedDict.items():
  20 			for orient, strata in slice.items():
  21 				for layer, point in strata.items():"""
  22 	from copy import deepcopy
  23 	import numpy as np
  24 	orient = {	'bone'  : np.zeros(3),
  25 				# 'muscle': np.zeros(3),
  26 				# 'fat'   : np.zeros(3),
  27 				'skin'	: np.zeros(3)	}
  28 	place = {	'Anterior' : deepcopy(orient),
  29 				'Posterior': deepcopy(orient),
  30 				'Medial'   : deepcopy(orient),
  31 				'Lateral'  : deepcopy(orient)	}
  32 	nestedDict={'Proximal' : deepcopy(place),
  33 				'Central'  : deepcopy(place),
  34 				'Distal'   : deepcopy(place)	}
  35 	return nestedDict
  36 
  37 
  38 def GetSubjectData(filename):
  39 	"""Reads in-vivo subject file MULTIS00X-1.xml
  40 		returns subjectData in a dictionary"""
  41 	subjectData = {}
  42 	tree = et.ElementTree(file=filename)
  43 	root = tree.getroot()
  44 	subject = root.find('Subject_Data')
  45 	measurements = subject.find('Anatomical_Measurements')
  46 	limbs = []
  47 	limbs.append(measurements.find('Lower_Leg') )
  48 	limbs.append(measurements.find('Lower_Arm'))
  49 	limbs.append(measurements.find('Upper_Leg'))
  50 	limbs.append(measurements.find('Upper_Arm'))
  51 	#
  52 	for limb in limbs:
  53 		subjectData[limb.tag]={}
  54 		for cluster in limb:
  55 			children = list(cluster)
  56 			descriptor, value = children[0].text, children[1].text
  57 			descriptor = descriptor.replace(' ','')
  58 			annoyingDescriptions =('Hip', 'Knee', 'Elbow', 'Shoulder')
  59 			for annoy in annoyingDescriptions:
  60 				descriptor=descriptor.replace(annoy,'')
  61 			subjectData[limb.tag][descriptor] = float(value)
  62 	return subjectData
  63 def GetCadaverData(filename):
  64 	"""Reads in-vitro cadaver subject file CMULTIS00X-1.xml
  65 		this file is structured slightly differently than in-vivo
  66 		returns subjectData in a dictionary"""
  67 	subjectData = {}
  68 	tree = et.ElementTree(file=filename)
  69 	root = tree.getroot()
  70 	subject = root.find('Subject_Data')
  71 	limbs = ['Arm', 'Leg']
  72 	segments = ['Upper', 'Lower']
  73 	measurements = []
  74 	for limb in limbs:
  75 		temp = subject.find('{}_-_Anatomical_Measurements'.format(limb))
  76 		for segment in segments:
  77 			measurements.append(temp.find('{}_{}'.format(segment,limb)))
  78 	for measure in measurements:
  79 		subjectData[measure.tag]={}
  80 		for cluster in measure:
  81 			children = list(cluster)
  82 			descriptor, value = children[0].text, children[1].text
  83 			descriptor = descriptor.replace(' ','')
  84 			annoyingDescriptions =('Hip', 'Knee', 'Elbow', 'Shoulder')
  85 			for annoy in annoyingDescriptions:
  86 				descriptor=descriptor.replace(annoy,'')
  87 			subjectData[measure.tag][descriptor] = float(value)
  88 	return subjectData
  89 
  90 def ReadUltrasoundThickness(limbSegment, xmlFile):
  91 	"""Read the xml to find thicknesses"""
  92 	import numpy as np
  93 	import os
  94 	import re
  95 	from lxml import etree as et
  96 	directory = os.path.dirname(xmlFile)
  97 	doc = et.parse(xmlFile)
  98 	root = doc.getroot()
  99 	bodyPart = root.find(limbSegment.replace('_', ''))
 100 	AllThicknesses = CreateNestedDictionary()
 101 	for loc in list(bodyPart):
 102 		place, orient = re.findall('[A-Z][a-z]*', loc.tag)  # Split on capitalization
 103 		thicknessFile = loc.get('Anatomical')
 104 		subDoc = et.parse(os.path.join(directory, thicknessFile))
 105 		subRoot = subDoc.getroot()
 106 		subject = subRoot.find("Subject").find("Source")
 107 		Thickness = {}
 108 		for frame in list(subject):
 109 			readings = frame.find('Thickness')
 110 			for value in list(readings):
 111 				try:				Thickness[value.tag.lower()].append(float(value.text))
 112 				except KeyError:	Thickness[value.tag.lower()] = [float(value.text)]
 113 		for t, v in Thickness.iteritems():
 114 			if np.all(np.isnan(v)):		# no data, use CT thickness evaluation instead
 115 				# except I don't have this for in-vivo data, only my test case...
 116 				pass
 117 				# v = []			# empty the list as it is only NaNs
 118 				# # note the follwing code is similar but NOT the same as the subDoc section above.
 119 				# subDoc = et.parse(os.path.join(directory, '../CTManual_zProbe',thicknessFile))	# adjust path to CT thickness directory
 120 				# subRoot = subDoc.getroot()
 121 				# subject = subRoot.find("Subject").find("Source")
 122 				# for frame in list(subject):
 123 				# 	value = frame.find('Thickness').find(t.capitalize())
 124 				# 	v.append(float(value.text))
 125 			elif  np.any(np.isnan(v)):		v = v[~np.isnan(v)]		# just use the good values
 126 			Thickness[t] = sum(v) / len(v)
 127 		AllThicknesses[place][orient] = Thickness
 128 	return AllThicknesses
 129 def ReadInVitroPositions(limbSegment, xmlFile):
 130 	"""Read the xml to find the in vitro ultrasound positions"""
 131 	import re
 132 	from lxml import etree as et
 133 	basePoints = CreateNestedDictionary()
 134 	baseDirections = CreateNestedDictionary()
 135 	newFile = xmlFile.replace('/','\\')
 136 	print newFile, os.path.isfile(newFile)
 137 	doc = et.parse(xmlFile)
 138 	root = doc.getroot()
 139 	bodyPart = root.find(limbSegment.replace('_',''))
 140 	for loc in list(bodyPart):
 141 		place = loc.tag
 142 		place, orient = re.findall('[A-Z][a-z]*', place)	# Split on capitalization
 143 		pos = loc.find("Anatomical").find("USPosition")
 144 		x = float( pos.find("x").get('value') ) * 1000
 145 		y = float( pos.find("y").get('value') ) * 1000
 146 		z = float( pos.find("z").get('value') ) * 1000
 147 		roll = float( pos.find("roll" ).get('value') )
 148 		pitch= float( pos.find("pitch").get('value') )
 149 		yaw  = float( pos.find("yaw"  ).get('value') )
 150 		basePoints[place][orient]['skin'] = [x, y, z]
 151 		baseDirections[place][orient] = findDirectionFromRPY(roll, pitch, yaw)
 152 	return basePoints, baseDirections
 153 
 154 def findDirectionFromRPY(roll, pitch, yaw, reverse=False):
 155 	from numpy import sin, cos, dot, array
 156 	r = roll
 157 	p = pitch
 158 	y = yaw
 159 	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)],
 160 								 [ cos(p)*sin(y),	sin(r)*sin(p)*sin(y) + cos(r)*cos(y),	cos(r)*sin(p)*sin(y) - sin(r)*cos(y)],
 161 								 [-sin(p), 			sin(r)*cos(p), 							cos(r)*cos(p)						]	])
 162 	if reverse:		dirVector = dot(transformMatrix, [0, 0,-1])
 163 	else:			dirVector = dot(transformMatrix, [0, 0, 1])  # ultrasound probe points in +Z dir
 164 	return dirVector
 165 
 166 def SimpleMorphing(renderer, subjectDir, cadaverDir, meshes, limb, activeMorphing=True, zPlane=False):
 167 	"""Morph a model (cadaver) into a shape matching an in-vivo subject"""
 168 	# draw meshes
 169 	myVTK.MakeStl(ren, meshes['bone'], color=colors.blanched_almond, opacity=0.1)
 170 	myVTK.MakeStl(ren, meshes['skin'], color=colors.cyan, opacity=0.1)
 171 	# Create handle for morphed meshes
 172 	morphs = {layer: mesh.replace('.','_morphed_{}.'.format(subjectDir.rsplit(os.sep,1)[1])) for layer,mesh in meshes.items() if layer is not 'bone'}
 173 	#
 174 	# Find files from positions relative to given directories
 175 	subjectFile =			os.path.join(subjectDir, 'Configuration', '{}.xml'.format(subjectDir.rsplit(os.sep, 1)[1]))
 176 	subjectUltrasound =		os.path.join(subjectDir, 'TissueThickness', 'UltrasoundManual', '{}_TA_inclusion.xml'.format(subjectDir.rsplit(os.sep, 1)[1]))
 177 	cadaverFile =			os.path.join(cadaverDir, 'Configuration', '{}.xml'.format(cadaverDir.rsplit(os.sep, 1)[1]))
 178 	cadaverUltrasound =		os.path.join(cadaverDir, 'TissueThickness', 'UltrasoundManual', '{}_TA_inclusion.xml'.format(cadaverDir.rsplit(os.sep, 1)[1]))
 179 	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())))
 180 	#
 181 	# Read cadaver subject file to get in-vitro experimentally meased values
 182 	cadaverMeasurements = GetCadaverData(cadaverFile)
 183 	cadaverLength  = cadaverMeasurements[limb]['Length'] * 10
 184 	cadaverCircum = {	'Central' : cadaverMeasurements[limb][ 'CentralCircumference']*10,
 185 						'Distal'  : cadaverMeasurements[limb][  'DistalCircumference']*10,
 186 						'Proximal': cadaverMeasurements[limb]['ProximalCircumference']*10	}
 187 	# Read subject file xml to find in-vivo experimentally measured sizes
 188 	subjectMeasurements = GetSubjectData(subjectFile)
 189 	subjectLength = subjectMeasurements[limb]['Length']*10
 190 	subjectCircum = {	'Central' : subjectMeasurements[limb][ 'CentralCircumference']*10,
 191 						'Distal'  : subjectMeasurements[limb][  'DistalCircumference']*10,
 192 						'Proximal': subjectMeasurements[limb]['ProximalCircumference']*10	}
 193 	#
 194 	# distance from central plane, file is distance from hip...so math is done now
 195 	subjectDistances= {	'Hip'     : subjectMeasurements[limb]['landmarktoCentralCircumference']*10,
 196 						'Proximal': subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectMeasurements[limb]['landmarktoProximalCircumference']*10,
 197 						'Central' : 0.0,
 198 						'Distal'  : subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectMeasurements[limb]['landmarktoDistalCircumference']*10,
 199 						'Knee'	  :	subjectMeasurements[limb]['landmarktoCentralCircumference']*10-subjectLength
 200 						}
 201 	#
 202 	# Read ultrasound thicknesses produces layer thickness
 203 	cadaverThicknesses = ReadUltrasoundThickness(limb, cadaverUltrasound)
 204 	subjectThickness = ReadUltrasoundThickness(limb, subjectUltrasound)
 205 	# sum layers for total thickness
 206 	cadaverThicknesses= {place:{orient:sum(strata.values()) for orient, strata in slice.items() } for place,slice in cadaverThicknesses.iteritems()}
 207 	subjectThickness  = {place:{orient:sum(strata.values()) for orient, strata in slice.items() } for place,slice in   subjectThickness.iteritems()}
 208 	#
 209 	# Read probe positions
 210 	probePoints, probeDirections = ReadInVitroPositions(limb, inVitroPositionsFile)
 211 	probeDirections = {place:{orient:strata for orient, strata in slice.items()} for place, slice in probeDirections.iteritems()}
 212 	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()}
 213 	#
 214 	if zPlane:	# if true forces probePoints into the same z plane
 215 		for place, slice in probePoints.iteritems():
 216 			for orient, strata in slice.items():
 217 				strata[2] = probePoints[place]['Anterior'][2]
 218 	#
 219 	# Find control point for morph
 220 	controlPoints = CreateNestedDictionary()
 221 	modelThicknesses = CreateNestedDictionary()
 222 	norms = {	'Anterior'	: np.array([ 0,-1, 0]),
 223 				'Lateral'	: np.array([-1, 0, 0]),
 224 				'Medial'	: np.array([ 1, 0, 0]),
 225 				'Posterior'	: np.array([ 0, 1, 0])	}
 226 	for place, slice in probePoints.iteritems():
 227 		# if place is 'Central':
 228 			for orient, strata in slice.items():
 229 				p1 = probePoints[place][orient]
 230 				p2 = p1 + norms[orient] * cadaverThicknesses[place][orient]
 231 				p2 = p1 + probeDirections[place][orient] * cadaverThicknesses[place][orient]
 232 				#
 233 				controlPoints[place][orient]['bone'] = myVTK.FindIntersectionOfStlAndLine(ren, meshes['bone'], p1, p2, colors.red  , visualize=True, embellish=True, extend=True)
 234 				controlPoints[place][orient]['skin'] = myVTK.FindIntersectionOfStlAndLine(ren, meshes['skin'], p1, p2, colors.green, visualize=False, embellish=True, extend=True)
 235 				modelThicknesses[place][orient] = myVTK.dist(controlPoints[place][orient]['bone'], controlPoints[place][orient]['skin'])
 236 	#
 237 	# myVTK.Show(ren, mouseHighlightActor=False)
 238 	#
 239 	morphedPoints = CreateNestedDictionary()
 240 	for place, slice in controlPoints.items():
 241 		for orient, strata in slice.items():
 242 			for layer, point in strata.items():
 243 				if layer is 'bone':		morphedPoints[place][orient][layer] = point
 244 				else:
 245 					if np.isnan(subjectThickness[place][orient]):	morphedPoints[place][orient][layer] = controlPoints[place][orient][layer]
 246 					else:	morphedPoints[place][orient][layer] = controlPoints[place][orient]['bone']-norms[orient]*subjectThickness[place][orient]
 247 					print place, orient, layer, '	',morphedPoints[place][orient][layer]
 248 	#
 249 	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'] )
 250 	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'] )
 251 	#
 252 	# Morph meshes with subprocess to speed things up.  Morphing of seperate layers will run concurrently.
 253 	if activeMorphing:		# This was used for testing purposes, it is left in because the user may want the output again, without recalculating the morph
 254 		print 'Morphing...',
 255 		processes = []
 256 		for layer, mesh in meshes.items():
 257 			if layer is not 'bone':
 258 				RunPygem.WritePygemParameterFile('{}.dat'.format(layer), controlArray,morphedArray)
 259 				processes.append( subprocess.Popen(['python', 'RunPygem.py', mesh, '{}.dat'.format(layer), morphs[layer]]) )
 260 		print 'this takes a while...',
 261 		for proc in processes:	proc.wait()
 262 		print 'finished.'
 263 	#
 264 	myVTK.MakeStl(ren, morphs['skin'], color=colors.orange, opacity=0.1)
 265 	#
 266 	morphThicknesses = CreateNestedDictionary()
 267 	for place, slice in morphedPoints.iteritems():
 268 		if place is 'Central':
 269 			for orient, strata in slice.items():
 270 				p1 = morphedPoints[place][orient]['skin']
 271 				p2 = p1 + norms[orient]*(modelThicknesses[place][orient]+subjectThickness[place][orient])/2
 272 				#
 273 				intersectBone = myVTK.FindIntersectionOfStlAndLine(ren, meshes['bone'], p1, p2, colors.red  , embellish=True, extend=True)	# no longer morphing bone
 274 				intersectSkin = myVTK.FindIntersectionOfStlAndLine(ren, morphs['skin'], p1, p2, colors.blue, embellish=True, extend=True, size=1)
 275 				morphThicknesses[place][orient] = myVTK.dist(intersectBone, intersectSkin)
 276 				diff = morphThicknesses[place][orient]-subjectThickness[place][orient]
 277 				
 278 	# myVTK.Show(ren, mouseHighlightActor=False)
 279 
 280 	#
 281 	planeNormal = (0,0,1)
 282 	planeSize = 250
 283 	inVitroColors = {'Proximal':colors.pink, 'Central':colors.green_pale, 'Distal':colors.cyan, 'Hip':colors.slate_grey, 'Knee':colors.yellow}
 284 	inVivo_Colors = {'Proximal':colors.red,  'Central':colors.green_dark, 'Distal':colors.blue, 'Hip':colors.black,		 'Knee':colors.orange}
 285 	#
 286 	modelCircum = {}
 287 	morphCircum = {}
 288 	for place, slice in controlPoints.items():
 289 		if place is 'Central':
 290 			cP = controlPoints[place]['Anterior']['skin']	# makes following lines somewhat more readable
 291 			mP = morphedPoints[place]['Anterior']['skin']
 292 			# myVTK.DrawPlane(ren, cP, planeNormal, size=planeSize, opacity=0.05, color=inVitroColors[place])  # for visualization only
 293 			plane = tvtk.Plane(origin=cP, normal=planeNormal)  # Create the plane as a function, that will be used to cut the stl.
 294 			modelCircum[place],c = myVTK.MakeStlCut(ren, meshes['skin'], plane, inVitroColors[place], opacity=0)
 295 			#
 296 			# myVTK.DrawPlane(ren, mP, normal, size=planeSize, opacity=0.05, color=inVivo_Colors[place])  # for visualization only
 297 			plane = tvtk.Plane(origin=mP, normal=planeNormal)  # Create the plane as a function, that will be used to cut the stl.
 298 			morphCircum[place],c = myVTK.MakeStlCut(ren, morphs['skin'], plane, inVivo_Colors[place], opacity=0)
 299 			# MBoneCircum[place],c = myVTK.MakeStlCut(ren, morphs['bone'], plane, colors.black, opacity=1)
 300 			# print '{:8s} Bone Morph circum {:6.2f}, Bone circum {:6.2f}'.format(place, MBoneCircum[place], bone_Circum[place])
 301 			print place
 302 			# print '	Model circum {: 6.2f}	thickness {: 6.2f}'.format(modelCircum[place], modelThicknesses[place])
 303 			print '	Morph	circum {: 6.2f}'.format(morphCircum[place]),
 304 			print '	Thickness:',
 305 			print ' Anterior {: 6.2f}'.format(morphThicknesses[place]['Anterior']),
 306 			print ' Posterior {: 6.2f}'.format(morphThicknesses[place]['Posterior']),
 307 			print ' Lateral {: 6.2f}'.format(morphThicknesses[place]['Lateral']),
 308 			print ' Medial {: 6.2f}'.format(morphThicknesses[place]['Medial'])
 309 			#
 310 			print '	Subj	circum {: 6.2f}'.format(subjectCircum[place]),
 311 			print '	Thickness:',
 312 			print ' Anterior {: 6.2f}'.format(subjectThickness[place]['Anterior'] ),
 313 			print ' Posterior {: 6.2f}'.format(subjectThickness[place]['Posterior']),
 314 			print ' Lateral {: 6.2f}'.format(subjectThickness[place]['Lateral']),
 315 			print ' Medial {: 6.2f}'.format(subjectThickness[place]['Medial'])
 316 			#
 317 			print '	Cadaver	circum {: 6.2f}'.format(cadaverCircum[place]),
 318 			print '	Thickness:',
 319 			print ' Anterior {: 6.2f}'.format(cadaverThicknesses[place]['Anterior'] ),
 320 			print ' Posterior {: 6.2f}'.format(cadaverThicknesses[place]['Posterior']),
 321 			print ' Lateral {: 6.2f}'.format(cadaverThicknesses[place]['Lateral']),
 322 			print ' Medial {: 6.2f}'.format(cadaverThicknesses[place]['Medial'])
 323 			#
 324 			print '	Model	circum {: 6.2f}'.format(modelCircum[place]),
 325 			print '	Thickness:',
 326 			print ' Anterior {: 6.2f}'.format(modelThicknesses[place]['Anterior'] ),
 327 			print ' Posterior {: 6.2f}'.format(modelThicknesses[place]['Posterior']),
 328 			print ' Lateral {: 6.2f}'.format(modelThicknesses[place]['Lateral']),
 329 			print ' Medial {: 6.2f}'.format(modelThicknesses[place]['Medial'])
 330 
 331 if __name__ == '__main__':
 332 	#
 333 	gender = 'female'
 334 	limb = 'Lower_Leg'	# choose from Upper_Arm, Upper_Leg, Lower_Arm, Lower_Leg
 335 	#
 336 	# This is a coded reminder of the model and invivo parings
 337 	if gender is 'female':
 338 		subjectDir = 'C:\Users\landisb\Workspace\InVivoSubjects\MULTIS037-1'
 339 		cadaverDir = 'C:\Users\landisb\Workspace\CadaverModels\CMULTIS008-1'
 340 	elif gender is 'male':
 341 		subjectDir = 'C:\Users\landisb\Workspace\InVivoSubjects\MULTIS094-1'
 342 		cadaverDir = 'C:\Users\landisb\Workspace\CadaverModels\CMULTIS006-1'
 343 	#
 344 	# 	Adjust meshes to match cadaver and limb segment choice
 345 	meshes = {
 346 		'bone': r'C:\Users\landisb\Downloads\Operation Multis-STL_Files\STL_Files\Lumped\008LL\meshed_bone_quad.stl',
 347 		'skin': r'C:\Users\landisb\Downloads\Operation Multis-STL_Files\STL_Files\Lumped\008LL\008LL_meshed_lump.stl'
 348 		}
 349 	#
 350 	ren = tvtk.Renderer(background=colors.white)
 351 	SimpleMorphing(ren, subjectDir, cadaverDir, meshes, limb, zPlane=False, activeMorphing=True)
 352 	print 'complete'
 353 	myVTK.Show(ren, mouseHighlightActor=False)