IsolatedAneurysmPP.py

import vtk
import numpy as np
from vtk.util import numpy_support as nps
import os
import codecs
import subprocess
import csv
import matplotlib.pyplot as plt


modelID = 'KDR13'
#SimModFlag = 'Rigid'
SimModFlag = 'FSI'
path = '/home/marsden-lab/Documents/alu2/' + modelID
AllResultsFileName = path +'/Results4MPA_720' + '/all_results.vtp'
AneurysmDirectory = path + '/Aneurysm/'
WallsMeshName = path + '/Mesh/mesh-complete/walls_combined.vtp'

#---------------------------------------------------------------------------------------------------------

def Get_TAWSS(Name, AneurysmDirectory, AllResultsFile, start, stop, incr):

	c=0

	numsteps = (stop - start)/incr + 1

	print 'Will be looking at ' + str(numsteps) + ' timesteps'
	print 'Reading walls ...'

	WallName =  Name.strip()
	FileName = AneurysmDirectory + WallName
	print FileName

	#Read vtp mesh files for the region of interest
	datareader = vtk.vtkXMLPolyDataReader()
	datareader.SetFileName(FileName)
	datareader.Update()
	data = datareader.GetOutput()

	#Get the number of points
	NoP = data.GetNumberOfPoints()
	print 'NoP in Get_TAWSS is ' + str(NoP)

	#Initialize total arrays
	total_tawss =np.zeros((1,NoP))
	total_osi =np.zeros((1,NoP))
	#total_point =np.zeros((NoP,3))
	#total_GlobalID = []

	#Read results file
	resultsreader = vtk.vtkXMLPolyDataReader()
	resultsreader.SetFileName(AllResultsFile)
	resultsreader.Update()
	results = resultsreader.GetOutput()

	#initialize arrays for tawss, osi, point (?)
	tawss_np = np.zeros((1,NoP))
	osi_np = np.zeros((1,NoP))
	time_wss = np.zeros((numsteps,NoP))
	final_wss = np.zeros(numsteps)
	point_np = np.zeros((NoP,3))


	#get the array of global node ids
	gNid_array = nps.vtk_to_numpy(results.GetPointData().GetArray('GlobalNodeID'))

	for k in range(NoP):
		gNid=data.GetPointData().GetArray('GlobalNodeID').GetTuple(k)
		#total_GlobalID.append( np.int(gNid[0]))

		#find indices in gnid_array where the node id matches our desired region
	   	index = np.where(gNid_array == np.int(gNid[0]))[0]
		tawss_np[0,k]=np.double(results.GetPointData().GetArray('vTAWSS').GetTuple(index[0]))
		osi_np[0,k]=np.double(results.GetPointData().GetArray('vOSI').GetTuple(index[0]))

		for t in range(numsteps):
			#print 'vWSS_0' + str(start + incr*t)
			time_wss[t,k] = np.linalg.norm(results.GetPointData().GetArray('vWSS_0'+str(start + incr*t)).GetTuple(index[0]))
			#temp_tawss[t] = np.double(results.GetPointData().GetArray('vWSS_0'+str(start + incr*t)).GetTuple(index[0]))

		#total_point[c,:] = np.double(data.GetPoints().GetPoint(k))
		total_tawss[0,c] = np.double(tawss_np[0,k])
		total_osi[0,c] = np.double(osi_np[0,k])
		#time_wss[0,c] = np.mean(temp_tawss)
		c=c+1

	for t in range(numsteps):
		final_wss[t] = np.mean(time_wss[t,:])

	tawss_vtk = nps.numpy_to_vtk(tawss_np[0])
	tawss_vtk.SetName('TAWSS')
	osi_vtk = nps.numpy_to_vtk(osi_np[0])
	osi_vtk.SetName('OSI')
	data.GetPointData().AddArray(tawss_vtk)
	data.GetPointData().AddArray(osi_vtk)

	data.GetCellData().RemoveArray('GlobalElementID2')
	datareader.Update()

	new=vtk.vtkXMLPolyDataWriter()
	new.SetInputData(data)
	new.SetFileName('Results_'+ WallName )
	new.Write()

	av_tawss = np.mean(total_tawss)
	av_osi = np.mean(total_osi)

	return (av_tawss, av_osi, final_wss)

#------------------------------------------------------------------------------------------------------


def Write_TAWSS(meshdata, total_tawss, total_GlobalID):

	NoMP = meshdata.GetNumberOfPoints()
	print 'in Write_TAWSS found NoMP = ' + str(NoMP)
	print 'total_GlobalID has len ' + str(len(total_GlobalID))
	meshTAWSS =np.zeros(NoMP)

	for k in range(NoMP):

		gMNid=meshdata.GetPointData().GetArray('GlobalNodeID').GetTuple(k)
		check = np.int(gMNid[0]) in total_GlobalID
		if check == True:
			index = total_GlobalID.index( np.int(gMNid[0]))
			meshTAWSS[k] = total_tawss[0,index]

	return meshTAWSS

#------------------------------------------------------------------------------------------------------
# function: WriteTStepData
# input:
#1) average tawss over entire region over time,
#2) average osi over entire region
#3) average tawss over region per time step

def WriteTStepData(av_tawss, av_osi, time_wss):
	t = 0
	for wss in time_wss:
		print str(t) +'     ' + str(wss).rjust(10)
		t = t+1
	plt.plot(range(20),time_wss)

#------------------------------------------------------------------------------------------------------

if __name__ == "__main__":

	#list all aneurysm mesh
	dummy = subprocess.Popen('ls ' + AneurysmDirectory,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
	dummylines = dummy.stdout.readlines()

	#get the total number of aneurysm
	NoC = len(dummylines)
	print 'NoC ' + str(NoC)

	#TODO logic to check that NoR + NoL = NoC? ? ?


	#read in the vtp wall combo mesh files
	meshreader = vtk.vtkXMLPolyDataReader()
	meshreader.SetFileName(WallsMeshName)
	meshreader.Update()
	meshdata = meshreader.GetOutput()

	for i in range(NoC):
		av_tawss, av_osi, time_wss= Get_TAWSS(dummylines[i], AneurysmDirectory, AllResultsFileName, 5000, 5950, 50)
		print '========================' + 'iteration #' + str(i) + '==============================='
		print 'aneurysm tawss (point average) = ' + str(av_tawss)
		print 'av_osi has value ' + str(av_osi)
		print '========================'
		print 'Preparing to give avg wss per time step: '
		WriteTStepData(av_tawss, av_osi, time_wss)