preprocessing.py
02915adb5f9bmaster

preprocessing.py
View Options

			"""
			This script preprocesses image and converts image to FE solver input file.
			"""
			from __future__ import print_function

			import sys
			import numpy as np

			import io_utils
			import utils
			from medtool import mic, fec


			def adjust_image_size(bone, config):
			"""
			Adjust image size to current FEelement size,
			that no layers of the image are removed Size in z direction has to fit.
			If empty layers are added in this dimension,
			this will create a weak layer at the bottom of the image.
			Expansions in x and y dimension,
			will probably not affect strength, but will lower stiffness.
			"""
			# get bone values
			BMD_array = bone["BMD_array"]
			MASK_array = bone["MASK_array"]
			Spacing = bone["Spacing"]

			# will get updated in subroutine
			FEelSize = np.copy(Spacing)
			# coarsening factor = FE element size / CT voxel size
			CoarseFactor = int(round(config["element_size"] / Spacing[0]))

			# Adjustment for BMD image and Mask
			BMD_array_adjusted = utils.adjust_image_size(
			BMD_array, CoarseFactor, utils.CropType.crop
			)
			# easier than changing subsquent code
			BMD_array = BMD_array_adjusted
			MASK_array_adjusted = utils.adjust_image_size(
			MASK_array, CoarseFactor, utils.CropType.crop
			)
			MASK_array = MASK_array_adjusted

			# set bone values
			bone["BMD_array"] = BMD_array
			bone["MASK_array"] = MASK_array
			bone["FEelSize"] = FEelSize
			bone["CoarseFactor"] = CoarseFactor

			return bone


			def calculate_bvtv(bone, IMTYPE):
			"""
			Calculate BVTV and mask images
			------------------------------
			Scaling, slope and intercept are printed out for review
			If image is already in BMD units, Scaling, Slope and
			Intercept are not applied. BVTVraw is scaled according to
			Hosseini et al. 2017
			(This scaling function could as well be defined externally).
			"""
			# get bone values
			Scaling = bone["Scaling"]
			Slope = bone["Slope"]
			Intercept = bone["Intercept"]
			BMD_array = bone["BMD_array"]

			print("\n ... prepare mask and BVTV images")
			print(" -> Scaling = ", Scaling)
			print(" -> Slope = ", Slope)
			print(" -> Intercept = ", Intercept)

			if IMTYPE.find("BMD") > -1:
			# if image is already in BMD units (e.g. Hosseinis data)
			BVTVraw = BMD_array / 1200.0
			elif IMTYPE.find("NATIVE") > -1:
			BMD_array = (BMD_array / Scaling) * Slope + Intercept
			BVTVraw = BMD_array / 1200.0 # if image is in native units

			# Scaling of BVTV 61um to BVTV 11.4um [Hosseini2017]
			# Seg_scaling_slope = 0.963
			# Seg_scaling_intercept = 0.03814

			BVTVscaled = BVTVraw
			# BVTVscaled = Seg_scaling_slope * BVTVraw
			# + Seg_scaling_intercept

			# set bone values
			bone["BVTVraw"] = BVTVraw
			bone["BVTVscaled"] = BVTVscaled

			return bone


			def generate_full_block_mesh(bone, filenames):
			# Medical Image Converter:
			# for converting segmented image to voxel mesh
			mymic = mic.mic()
			# Finite Element Converter:
			# for reading and writing Abaqus .inp and Ensight .case files
			myfec = fec.fec()
			# Medical Image Analysis:
			# for the fabric calculations
			# mymia = mia.mia()

			# get bone values
			BVTVscaled = bone["BVTVscaled"]
			MASK_array = bone["MASK_array"]
			FEelSize = bone["FEelSize"]
			CoarseFactor = bone["CoarseFactor"]

			BVTVb_masked = np.copy(BVTVscaled)
			BVTVb_masked[MASK_array == 0] = 0

			# The original mesh is a full block from which
			# the elements without bone are going to be removed
			MESH = np.ones_like(BVTVscaled)
			# Beware, Spacing gets updated...
			MESH = mymic.changeResolutionF(MESH, CoarseFactor, FEelSize)
			# weird np array convention (z,y,x)
			MESH = MESH.transpose(2, 1, 0)
			# esize = Spacing # real element size in mm
			# *****************************************************************

			io_utils.print_mem_usage()

			# *****************************************************************
			# Generate full block mesh
			# ----------------------------------------------------------------
			print("\n ... generate full block mesh (Abaqus inp file)")
			INPname = filenames["INPname"]

			mymic.writeAbaqusGeneral(INPname, MESH, FEelSize, templateFile=None, smooth=None)
			inp = myfec.readAbaqus(INPname)
			# title = inp[0]
			nodes = inp[1]
			# nsets = inp[2]
			elems = inp[3]
			elsets = inp[4]
			# element set "SET1" correspondes to the full block
			# with the size of the image
			# A new element set "BONE" is created after material mapping
			elems = {elem: elems[elem] for elem in elsets["SET1"]}
			nodes = {node: nodes[node] for elem in elems for node in elems[elem].get_nodes()}
			elsets = {}
			print(" - finished")
			# *****************************************************************

			# set bone values
			bone["elems"] = elems
			bone["nodes"] = nodes
			bone["elsets"] = elsets
			bone["MESH"] = MESH

			return bone


			def calculate_fabric(bone, config):
			"""
			Compute isotropic fabric
			"""

			FTYPE = config["ftype"]
			if FTYPE.find("Iso") > -1 or FTYPE.find("iso") > -1 or FTYPE.find("test") > -1:
			print("\n ... compute isotropic fabric")
			evalue, evect = utils.compute_isoFAB()
			print(" - finished")

			io_utils.print_mem_usage()

			bone["evalue"] = evalue
			bone["evect"] = evect

			return bone


			def material_mapping(bone, config, umat_parameters, filenames):
			"""
			Material Mapping
			"""

			# get bone values
			BVTVscaled = bone["BVTVscaled"]
			MASK_array = bone["MASK_array"]
			FEelSize = bone["FEelSize"]
			Spacing = bone["Spacing"]
			elems = bone["elems"]
			nodes = bone["nodes"]
			elsets = bone["elsets"]
			evalue = bone["evalue"]
			evect = bone["evect"]

			# TODO change for new phase BVTV
			elsets["BONE"] = []
			RHOb = {}
			RHOb_FE = {}
			RHOb_corrected = {}
			PHIb = {}
			mm = {}
			m = {}
			cogs = {}
			superp_stiff = {}
			superp_comp = {}
			superp_ffff = {}
			superp_ff = {}

			ROI_BVTV_size = config["roi_bvtv_size"]

			for i, elem in enumerate(elems):
			cog = np.mean(
			[np.asarray(nodes[node].get_coord()) for node in elems[elem].get_nodes()],
			axis=0,
			)

			PHIbone, Xb, Yb, Zb = utils.computePHI(cog, Spacing, FEelSize[0], MASK_array)

			if PHIbone > 0.0:
			elsets["BONE"].append(elem)
			BVTVbone = utils.computeBVTV_onephase(
			cog, Spacing, ROI_BVTV_size, BVTVscaled, MASK_array, PHIbone
			)
			BVTVbone_FE = utils.computeBVTV_FEel(
			cog, Spacing, FEelSize, BVTVscaled, MASK_array
			)

			RHOb[elem] = BVTVbone
			RHOb_FE[elem] = BVTVbone_FE
			PHIb[elem] = PHIbone
			RHOb_corrected[elem] = BVTVbone * PHIbone

			m[elem], mm[elem] = evalue, evect
			cogs[elem] = cog # included just for debugging reasons

			# material superposition
			# ----------------------
			factor = 0.05
			if BVTVbone <= factor:
			BVTVbone = factor
			SSSSt, CCCCt, FFFFt, FFt = utils.trab_ssss(
			umat_parameters, PHIbone, BVTVbone, evalue[0], evalue[1], evalue[2]
			)
			superp_stiff[elem] = SSSSt
			superp_comp[elem] = CCCCt
			superp_ffff[elem] = FFFFt
			superp_ff[elem] = FFt

			sys.stdout.write(
			"\r" + " ... material mapping element " + str(i) + "/" + str(len(elems))
			)
			sys.stdout.flush()

			# conversion to np array for calculation of BVTV
			RHOb_array = np.array(RHOb.values())
			# RHOb_corrected_array = np.array(RHOb_corrected.values())
			RHOb_FE_array = np.array(RHOb_FE.values())

			PHIb_array = np.array(PHIb.values())

			# superp_stiff_array = np.array(superp_stiff.values())
			# superp_comp_array = np.array(superp_comp.values())
			# superp_ffff_array = np.array(superp_ffff.values())
			# superp_ff_array = np.array(superp_ff.values())

			elems = {elem: elems[elem] for elem in elsets["BONE"]}
			nodes = {node: nodes[node] for elem in elems for node in elems[elem].get_nodes()}

			myfec = fec.fec()
			INPname = filenames["INPname"]
			myfec.writeAbaqus(INPname, None, nodes, None, elems, elsets, NscaResults=None)
			# *****************************************************************
			marray = np.real(np.mean([np.asarray(m[elem]) for elem in m.keys()], axis=0))
			mmarray1 = np.real(np.mean([np.asarray(mm[elem][0]) for elem in m.keys()], axis=0))
			mmarray2 = np.real(np.mean([np.asarray(mm[elem][1]) for elem in m.keys()], axis=0))
			mmarray3 = np.real(np.mean([np.asarray(mm[elem][2]) for elem in m.keys()], axis=0))

			io_utils.print_mem_usage()

			print("... material mapping finished")

			# set bone values
			bone["PHIb"] = PHIb
			bone["RHOb"] = RHOb
			bone["RHOb_array"] = RHOb_array
			bone["RHOb_FE_array"] = RHOb_FE_array
			bone["PHIb_array"] = PHIb_array
			bone["elems"] = elems
			bone["nodes"] = nodes
			bone["elsets"] = elsets
			bone["cogs"] = cogs
			bone["marray"] = marray
			bone["mmarray1"] = mmarray1
			bone["mmarray2"] = mmarray2
			bone["mmarray3"] = mmarray3
			bone["m"] = m
			bone["mm"] = mm
			bone["superp_stiff"] = superp_stiff
			bone["superp_comp"] = superp_comp
			bone["superp_ffff"] = superp_ffff
			bone["superp_ff"] = superp_ff

			return bone


			def set_summary_variables(bone):
			"""
			Computes variables for summary file
			"""

			# get bone values
			BMD_array = bone["BMD_array"]
			BVTVscaled = bone["BVTVscaled"]
			MASK_array = bone["MASK_array"]
			FEelSize = bone["FEelSize"]
			RHOb_array = bone["RHOb_array"]
			RHOb_FE_array = bone["RHOb_FE_array"]
			PHIb_array = bone["PHIb_array"]

			variables = {}

			# For entry in summary file
			# variables['BVTVscaling_string'] = "BVTV scaled = {:.2f}
			# * BVTVraw + {:.2f}".format(
			# Seg_scaling_slope, Seg_scaling_intercept)

			# Computation of variables for summary file
			variables["apparentBMD"] = np.mean(BMD_array[MASK_array != 0])
			variables["Full_Mask_Volume"] = np.sum(MASK_array > 0) * FEelSize[1] ** 3
			variables["apparentBVTV"] = variables["apparentBMD"] / 1200
			variables["BV_image"] = np.mean(BVTVscaled[MASK_array != 0]) * float(
			np.count_nonzero(MASK_array) * FEelSize[1] ** 3
			)

			print(
			"\n Apparent BMD of grayscale image: {:.2f} mgHA/ccm".format(
			variables["apparentBMD"]
			)
			)

			# mean tissue BVTV
			variables["ROI_tissue_BVTV"] = np.sum(RHOb_array * PHIb_array) / np.sum(PHIb_array)

			# mean simulation BVTV
			variables["ROI_simulation_BVTV"] = np.sum(RHOb_array * PHIb_array) / len(RHOb_array)

			# mean tissue BVTV FE elements
			variables["FE_tissue_BVTV"] = np.sum(RHOb_FE_array * PHIb_array) / np.sum(
			PHIb_array
			)

			# mean simulation BVTV FE
			variables["FE_simulation_BVTV"] = np.sum(RHOb_FE_array * PHIb_array) / len(
			RHOb_array
			)

			# BMC summarized from ROI
			variables["ROI_BMC_elements"] = (
			np.sum(RHOb_array * PHIb_array * FEelSize[1] ** 3) * 1200 * 0.001
			)

			# BMC summarized from FE
			variables["FE_BMC_elements"] = (
			np.sum(RHOb_FE_array * PHIb_array * FEelSize[1] ** 3) * 1200 * 0.001
			)

			variables["BMC_ROI_to_FE"] = (
			variables["ROI_BMC_elements"] / variables["FE_BMC_elements"]
			)

			# Mask volume
			variables["ROI_Mask_Volume"] = np.sum(PHIb_array * FEelSize[1] ** 3)

			# Parameters calculated on global image
			variables["bone_mass_mgHA"] = (
			variables["apparentBMD"] * variables["Full_Mask_Volume"] * 0.001
			)
			variables["BMC_image_to_FE"] = (
			variables["bone_mass_mgHA"] / variables["FE_BMC_elements"]
			)
			variables["volume_ROI_to_FE"] = (
			variables["ROI_Mask_Volume"] / variables["Full_Mask_Volume"]
			)

			io_utils.print_mem_usage()

			return variables

bonepy/bonepy/preprocessing.py02915adb5f9bmaster

preprocessing.pyView Options

bonepy/bonepy/preprocessing.py
02915adb5f9bmaster

preprocessing.py
View Options