-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA256

Hello,

I'm attempting to replicate the study from Figure 14 here: https://doi.org/10.1016/j.media.2019.06.004 with a regular mesh, not using the meshless method.

I have:

• A masked brain that can be used to generate a mesh from voxels
• The location of points to set for the displacement relative to points on the brain mesh
• Coordinates and triangles connecting the coordinates for the inner skull surface

I would like to get out of the model:

• The displaced vertices so that the original brain image can be resampled

What I'm stuck on:

• How to define the boundary conditions per coordinate so that the correct coordinates are displaced
• How to get the displacement vectors out of the solved problem

Bonus:

• If there is a way to do this with the inner skull surface acting as a fixed boundary condition frictionless with the brain, that would be preferable, I'm a more lost on how to do this though

If someone had time to point me in the right direction a bit, that would be greatly appreciated. I think if I knew a few more of the correct functions to use, I would be very close. Below is what I have so far (requires pip install mne).

import os.path as op
import numpy as np

import mne
import nibabel as nib

from sfepy.base.base import IndexedStruct
from sfepy.discrete import (FieldVariable, Material, Integral,
                            Equation, Equations, Problem)
from sfepy.discrete.fem import FEDomain, Field
from sfepy.mesh.mesh_generators import gen_mesh_from_voxels
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.postprocess.viewer import Viewer
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson

misc_path = mne.datasets.misc.data_path()

# get original and projected positions
raw = mne.io.read_raw(op.join(misc_path, 'ecog', 'sample_ecog_ieeg.fif'))
trans = mne.coreg.estimate_head_mri_t(
    'sample_ecog', op.join(misc_path, 'ecog'))
montage = raw.get_montage()
montage.apply_trans(trans)  # transform to MRI coordinates
ch_pos_disp = np.array([
    pos for pos in montage.get_positions()['ch_pos'].values()])

raw.info = mne.preprocessing.ieeg.project_sensors_onto_brain(
    raw.info, trans, 'sample_ecog', subjects_dir=op.join(misc_path, 'ecog'))
montage = raw.get_montage()
montage.apply_trans(trans)  # transform to MRI coordinates
ch_pos_orig = np.array([
    pos for pos in montage.get_positions()['ch_pos'].values()])

# define model
dist = 10  # displace all vertices up to 1 cm from points
brainmask = nib.load(op.join(misc_path, 'ecog', 'sample_ecog',
                             'mri', 'brainmask.mgz'))

mesh = gen_mesh_from_voxels(brainmask, [1, 1, 1])
domain = FEDomain('domain', mesh)

omega = domain.create_region('Omega', 'all')
fixed = domain.create_region(
    'Fixed', XXX, 'facet')
disp = domain.create_region(
    'Displaced', XXX, 'facet')

field = Field.from_args('fu', np.float64, 'vector', omega,
                        approx_order=2)

u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')

brain = Material('Brain', D=stiffness_from_youngpoisson(
    dim=3, young=3000, poisson=0.49), density=1000)
# csf = Material('Brain', D=stiffness_from_youngpoisson(
#    dim=3, young=100, poisson=0.49), density=1000)

integral = Integral('i', order=1)

term = Term.new('dw_lin_elastic(m.D, v, u)',
                integral, omega, m=brain, v=v, u=u)
eq = Equation('balance', term)
eqs = Equations([eq])

fix_u = EssentialBC('fix_u', fixed, {'u.all': 0.0})
disp_u = EssentialBC(
    'disp_u', disp, {'u.0': XXX, 'u.1': XXX, 'u.2': XXX})

ls = ScipyDirect({})

nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)

pb = Problem('deform', equations=eqs)

# view setup
'''
pb.save_regions_as_groups('regions')
view = Viewer('regions.vtk')
view()
'''

pb.set_bcs(ebcs=Conditions([fix_u, disp_u]))

pb.set_solver(nls)

status = IndexedStruct()
state = pb.solve(status=status)

print('Nonlinear solver status:\n', nls_status)
print('Stationary solver status:\n', status)

pb.save_state('brain_deform.vtk', state)

view = Viewer('brain_deform.vtk')
view(vector_mode='warp_norm', rel_scaling=1,
     is_scalar_bar=True, is_wireframe=True)

# modify original image
T1 = nib.load(op.join(misc_path, 'ecog', 'sample_ecog', 'mri', 'T1.mgz'))
# TO DO: get displacement vector from each vertex out, resample image
# based on T1 image intensity of original element positions

Thanks,

Alex Rockhill

-----BEGIN PGP SIGNATURE----- Version: BCPG v1.65

iQEcBAABCgAGBQJhjq3nAAoJEOCF7T6SI2rbi08H/j9PIYQxCaT0aDe5+EmAKmaw I+rZTM+gKBoDZkJD+XTUqUdXKALLxqdpNShi5JAy5MTeXS+CGJ9PufoL3DgDoWJE czfbczqGXzZStsKL7qVROEBF7sP2uvWNxdVYTXqUWVxXdn814SBn+pGHxPF7JShU XzTP1D+Bp/mlPL7/3rcw/qpI0LnFBdLKIcnwEERetEfpjKdj1v0/sAOpAJknHhGO yKwU+4M5bt95Q29S1qyzM5ybMFToj8SEdEBa8qPgtdBt67NlYv4bfrpsgt79i/eq +xgbcWV5M6zDeJyhtlm7YTRhc9cCsxTdug0jB5TdDXWFFsy2kS1NNtacnsy7BVE= =Ntng -----END PGP SIGNATURE-----