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


276
views
2
6 months ago by
Ryan W  
I am trying to implement a simple advection-diffusion equation, similar to what is shown in ft08_navier_stokes_cylinder.py. I have 3D images of a patient's brain, which let's me measure the intial conditions (contrast agent) and the material properties (effective diffusivity), as a 3D matrix. The goal of this work is to use the imaging data to constrain the problem, and iteratively run the FEM with differing scaling parameters. I would like to be assign the ICs and material properties as an interpolation of the imaging measurements. The image analysis and mesh assembly is done in MATLAB, so I can simply read in the images.

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

Community: FEniCS Project

1 Answer


4
6 months ago by
Ryan W  
I have figured it out. Instead of doing the above, I simply inherit from the Expression class, and use scipy's RegularGridInterpolater.

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​
Hi Ryan,
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,
written 8 days ago by Ricardo Tenorio  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »