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
Hi Alex,
On 12/11/2021 19:09, Alexander Rockhill via SfePy wrote:
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
So you want to fix the brain points that are given by ch_pos_orig
and
displace the points, using the boundary conditions, given by ch_pos_disp
in
your script?
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
By this you mean having the inner brain surface sliding along a given surface (skull) using some kinematic constraints? That is currently not possible.
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
).
FWIW, I had to do pip install mne pooch nibabel tqdm
. Then I managed to run
the script, up to the mesh creation, after changing:
mesh = gen_mesh_from_voxels(brainmask, [1, 1, 1])
to
mesh = gen_mesh_from_voxels(brainmask.dataobj, [1, 1, 1])
You might also be interested in the sfepy.mesh.mesh_tools.smooth_mesh()
function.
r.
Exactly, I would like to fix the brain points that are given by ch_pos_orig
and displace the points, using the boundary conditions, given by ch_pos_disp
. Then, it would be great to get the displaced vertices out.
Thanks for the response and sorry that the dependancies weren't all installed by pip, thanks for the heads up.
Alex
On 11/12/2021 2:07 PM Robert Cimrman <cimrman3@ntc.zcu.cz> wrote:
Hi Alex,
On 12/11/2021 19:09, Alexander Rockhill via SfePy wrote:
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
So you want to fix the brain points that are given by
ch_pos_orig
and displace the points, using the boundary conditions, given bych_pos_disp
in your script?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
By this you mean having the inner brain surface sliding along a given surface (skull) using some kinematic constraints? That is currently not possible.
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
).FWIW, I had to do
pip install mne pooch nibabel tqdm
. Then I managed to run the script, up to the mesh creation, after changing:mesh = gen_mesh_from_voxels(brainmask, [1, 1, 1])
to
mesh = gen_mesh_from_voxels(brainmask.dataobj, [1, 1, 1])
You might also be interested in the
sfepy.mesh.mesh_tools.smooth_mesh()
function.r.
SfePy mailing list  sfepy@python.org To unsubscribe send an email to sfepyleave@python.org https://mail.python.org/mailman3/lists/sfepy.python.org/ Member address: aprockhill@mailbox.org
Alex Rockhill
Hi Alex,
On 13/11/2021 00:41, Alexander Rockhill via SfePy wrote:
Exactly, I would like to fix the brain points that are given by
ch_pos_orig
and displace the points, using the boundary conditions, given bych_pos_disp
. Then, it would be great to get the displaced vertices out.
To define the regions, use something like:
from sfepy.discrete import Function, Functions from scipy.spatial import cKDTree
def _get_fixed(coors, domain=None): print(ch_pos_orig) kdtree = cKDTree(coors) cc = kdtree.query(ch_pos_orig) ... flag = ... # True where coors in the region. return np.nonzero(flag)[0]
get_fixed = Function('get_fixed', _get_fixed) fixed = domain.create_region( 'Fixed', 'vertices by get_fixed', 'facet', functions=Functions([get_fixed]), )
I could not finish it, because the points are scaled differently than the mesh coordinates, but you get the idea: the kdtree query should return the vertices in the mesh closest to the points. Analogously for the displaced region.
Concerning the boundary conditions, I am not sure what you want  also describe the displacements using a function of coordinates?
Exporting the displacements should be easy, just reshape the state:
u = state().reshape((1, 3))
Then you can index u with a region's vertices.
Thanks for the response and sorry that the dependancies weren't all installed by pip, thanks for the heads up.
No problem, pip is bad at resolving dependencies.
r.
Alex
On 11/12/2021 2:07 PM Robert Cimrman <cimrman3@ntc.zcu.cz> wrote:
Hi Alex,
On 12/11/2021 19:09, Alexander Rockhill via SfePy wrote:
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
So you want to fix the brain points that are given by
ch_pos_orig
and displace the points, using the boundary conditions, given bych_pos_disp
in your script?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
By this you mean having the inner brain surface sliding along a given surface (skull) using some kinematic constraints? That is currently not possible.
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
).FWIW, I had to do
pip install mne pooch nibabel tqdm
. Then I managed to run the script, up to the mesh creation, after changing:mesh = gen_mesh_from_voxels(brainmask, [1, 1, 1])
to
mesh = gen_mesh_from_voxels(brainmask.dataobj, [1, 1, 1])
You might also be interested in the
sfepy.mesh.mesh_tools.smooth_mesh()
function.
participants (2)

Alexander Rockhill

Robert Cimrman