### Is there a way to assign an interpolation of an image as initial conditions or material properties?

My mesh is 3D tetrahedral, in the shape of a patient's brain. The mesh is 2x higher resolution than the IC and material property images I have. Currently, I assign material properties based on proximity to the nearest voxel (3D pixel) center. I feel like this is brute-forcing the problem, and that there is a more elegant and efficient way to simply assign the measured images as an interpolated field.

My current strategy so far:

```
from fenics import *
import time
import scipy.io as spio #library for reading MATLAB .MAT files
import numpy as np
from tqdm import tqdm #console progress bar, with timing estimates
#define constants
T= 20.0
num_steps = 20
dt= T / num_steps
tol=0.1
class MeshProps:
#Parse .MAT storage file which stores mesh properties
def __init__(self,MATLAB):
raw_data=spio.loadmat(MATLAB)
props=raw_data['MeshProps']
self.ADC=props['ADC'][0,0] #Effective Diffusion Coeffient. Shape --> [N_voxels]
self.IC=props['IC'][0,0] #initial conditions. Shape --> [N_voxels]
self.VoxelCenters=props['VoxelCenters'][0,0] #center of voxels, for assinging
#above properties. Shape --> [3,N_voxels]
#Resolution is 1 unit per pixel, from center to center
#mesh is defined from this spatial mapping
class Voxel(SubDomain):
#Create a subdomain for each individual voxel
#Assign each element the material properties for nearest voxel
def __init__(self, VoxelCenters, CurrentIdx):
self.VoxelCenters=VoxelCenters
self.CurrentIdx=CurrentIdx
SubDomain.__init__(self)
def inside(self, x, on_boundary):
#Label points inside a voxel if they are inside the voxel
flagx = np.abs(x[0]-self.VoxelCenters[self.CurrentIdx,0]) <= 0.5 + tol
flagy = np.abs(x[1]-self.VoxelCenters[self.CurrentIdx,1]) <= 0.5 + tol
flagz = np.abs(x[2]-self.VoxelCenters[self.CurrentIdx,2]) <= 0.5 + tol
return np.logical_and(np.logical_and(flagx,flagy),flagz)
class VoxelWiseParameter(Expression):
def __init__(self,materials,param_array,**kwargs):
self.materials=materials
self.param_array=param_array
def eval_cell(self,values,x,cell):
values[0]=self.param_array[self.materials[cell.index]]
#Load in the mesh
mesh=Mesh('patientbrain_mesh.xml')
mesh=refine(mesh)
#define the materials space
materials=CellFunction('size_t',mesh)
#Use the above defined class to parse MATLAB data
meshprops=MeshProps('params_byhand.mat')
for i in tqdm(range(0,max(meshprops.VoxelCenters.shape)),desc='Assigning BCs'):
voxel=Voxel(meshprops.VoxelCenters,i) #add a new subdomain for each voxel
voxel.mark(materials, i) #mark the elements with its indicator
#Save the materials so we don't have to run previous loop again
File('pat11_materials.xml.gz') << materials
#assign ICs
u_0=VoxelWiseParameter(materials,meshprops.IC,degree=0)
#Assign material properties
adc=VoxelWiseParameter(materials, meshprops.ADC,degree=0) #apparent diffusion constant
D=Constant(0.001) #base diffusion, to be varied iteratevly (eventually)
#Set up the FEM
V=FunctionSpace(mesh, 'Lagrange', 1)
u_n = interpolate(u_0, V)
u_n.rename('u', 'initial value')
vtkfile = File('ModelOutput.pvd')
vtkfile << (u_n,0.0)
u=TrialFunction(V)
v=TestFunction(V)
F=u*v*dx + D*adc*dt*dot(grad(u),grad(v))*dx
a,L = lhs(F), rhs(F)
A=assemble(a)
u=Function(V)
u.rename('u','solution')
t=0
#Implement the Transient Model
for n in tqdm(range(num_steps),desc='Transient Model'):
t += dt
B=assemble(L)
solve(A, u.vector(), B, 'bicgstab', 'hypre_amg')
vtkfile << (u,float(t))
u_n.assign(u)
#Display a warning if a large portions of material remain unmarked
if np.count_nonzero(materials.array()==0) > 50:
print('Warning: Not all all elements were assigned a material.')
```

Using the above method, I often trigger the end warning. After visual inspection of the materials, it looks like corners, edges, and regions near the edge of the mesh remain unmarked. This code throws no errors, other than the one I hard-coded in. The initial conditions, when viewed in Paraview, remain sparse, with holes in the mesh missing. However, the simulation yeilds results based on the differing material properties assigned.

I've attached an image of the mesh, as well as the final result of the mesh implementation.

TD;DR Is there a way to assign an interpolation of an image as initial conditions or material properties? I have brute-forced it, and apply a different material to each individual imaging voxel. It seems to not throw errors, but it leaves significant portions of the mesh unmarked.

Many thanks ahead of time!

Ryan

### 1 Answer

Applying the interpolation still takes a while (~30 minutes), so I may look into a way to do this in C++. Otherwise, this appears to work.

```
from scipy.interpolate import RegularGridInterpolator
class InterpolatedParameter(Expression):
def __init__(self,X,Y,Z,param_image,**kwargs):
self.X=X #a numpy array giving the Z-spacing of the image
self.Y=Y #same for Y
self.Z=Z #same for Z
self.param_image=param_image #The image of measured material property
def eval_cell(self,values,x,cell):
interp_handle=RegularGridInterpolator((self.X,self.Y,self.Z),self.param_image)
values[0]=interp_handle(x)
adc=InterpolatedPatameter(range(0,param_image.shape[0]),\
range(0,param_image.shape[1]),\
range(0,param_image.shape[2]),\
param_image,\ #a 3d image of aparent diffusivity
degree=0)
adc=interpolate(adc,V) #use the fenics interpolate function to apply the interpolation to the DOFs
```

could you upload a sample of param_image? I need to use a 2D array as initial condition and I think your answer is the closest way to do what I am looking for,

thanks a lot,