Coverage for / builds / tms-localization / papers / tmsloc_proto / scripts / 02_make_msh_from_mri.py: 63%
449 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-01 07:48 +0000
1import os
2import sys
3import time
4import h5py
5import shutil
6import pynibs
7import simnibs
8import warnings
9import argparse
10import datetime
11import numpy as np
12import nibabel as nib
14from simnibs import __version__ as sim_v
15from pynibs.util.simnibs_io import smooth_mesh
17# Get command line arguments
18parser = argparse.ArgumentParser(description='Creates a FEM mesh from MRI data')
19parser.add_argument('-s', '--fn_subject', help='Path to *.hdf5-file of subject', required=True)
20parser.add_argument('--charm_ini', help='charm.ini file to use.',
21 default=f"{os.path.dirname(os.path.realpath(__file__))}/configs/charm_410.ini")
22parser.add_argument('-m', '--mesh_idx', help='The mesh name from subject object', required=True)
23parser.add_argument('-l', '--layer', help='Create cortical layer', action='store_true', required=False,
24 default=False)
25parser.add_argument('--charm_cmd_params', help='Additional parameters to charm, e.g. forcerun or mesh',
26 required=False,
27 default='')
28parser.add_argument('--mesh', help='Skip segmentation and only mesh headmodel', action='store_true',
29 required=False,
30 default=False)
31parser.add_argument('--debug', help='Run charm in debug mode saving additional intermediate output',
32 action='store_true', required=False, default=False)
34args = parser.parse_args()
35charm_cmdline = "charm"
36dwi2cond_cmdline = "dwi2cond"
37if os.system(f"which {dwi2cond_cmdline}") == 256:
38 simnibs_dir = os.path.split(simnibs.__file__)[0]
39 dwi2cond_cmdline = f"{simnibs_dir}/external/dwi2cond"
40 assert os.system(f"which {dwi2cond_cmdline}") == 0, f"Cannot find dwi2cond executable."
41meshmesh_cmdline = "meshmesh"
42layer = args.layer
43mesh_only = args.mesh
45debug_flag = ""
46if args.debug:
47 debug_flag = "--debug"
50def logprint(text):
51 """Just a wrapper to print and write text to a logfile."""
52 with open(log_file, "a") as logfile:
53 logfile.write(f'{text}\n')
54 print(text)
57fn_subject = os.path.abspath(args.fn_subject)
58subject = pynibs.load_subject(fn_subject)
59mesh_idx = args.mesh_idx
60charm_ini_fn = args.charm_ini
62assert mesh_idx in subject.mesh, f"Mesh '{mesh_idx}' not found in {args.fn_subject}"
63mesh_folder = subject.mesh[mesh_idx]["mesh_folder"]
64log_file = os.path.join(mesh_folder, f"info_{datetime.datetime.now():%Y%m%d_%H%M%S}.log")
65charm_cmd_params = args.charm_cmd_params
66if charm_cmd_params != '' and not charm_cmd_params.startswith('--'):
67 charm_cmd_params = '--' + charm_cmd_params
69try:
70 use_fs = subject.mesh[mesh_idx]["use_fs"]
71except KeyError:
72 use_fs = False
75def print_mesh_info(fn):
76 logprint("- " * 32)
77 logprint(f"Mesh information: {fn}")
78 pref = os.path.basename(fn)
79 with h5py.File(fn, 'r') as f: # f = h5py.File(mesh_fn_smoothed_hdf5,'r')
80 logprint(f"{pref}: total nodes: {f['/mesh/nodes/node_number'].shape[0]: >7}")
81 logprint(f"{pref}: total elms: {f['/mesh/elm/elm_number'].shape[0]: >7}")
82 logprint(f"{pref}: total tris: {np.sum(f['/mesh/elm/elm_type'][:] == 2): >7}")
83 logprint(f"{pref}: total tets: {np.sum(f['/mesh/elm/elm_type'][:] == 4): >7}")
84 logprint("- " * 32)
85 for t in np.unique(f['/mesh/elm/tet_tissue_type']):
86 logprint(f"{pref}: tet tissue type {t:0>4}: {np.sum(f['/mesh/elm/tet_tissue_type'] == t): >7}")
87 logprint("- " * 32)
88 for t in np.unique(f['/mesh/elm/tri_tissue_type']):
89 logprint(f"{pref}: tri tissue type {t:0>4}: {np.sum(f['/mesh/elm/tri_tissue_type'] == t): >7}")
90 logprint("- " * 32)
93time_mesh = 0
94time_roi = 0
95time_dwi = 0
97# logprint simulation info
98# ================================================================================
99os.makedirs(os.path.split(log_file)[0], exist_ok=True)
100logprint("-" * 64)
101logprint(f"{parser.description}")
102args_dict = vars(args)
103for key in args_dict.keys():
104 logprint(f"{key: >15}: {args_dict[key]}")
105logprint("-" * 64)
107if charm_ini_fn:
108 # this will always be appended to all charm, meshmesh calls
109 charm_ini_fn = f" --usesettings={os.path.abspath(charm_ini_fn)}"
110else:
111 charm_ini_fn = ''
113# load subject info
114mesh_dict = subject.mesh[mesh_idx]
115mri_idx = mesh_dict["mri_idx"]
117# set approach specific options
118approach = mesh_dict["approach"]
120if approach == "charm":
121 vertex_density = None
122 numvertices = None
123else:
124 raise NotImplementedError(f"Approach {approach} not implemented. Use 'charm' or switch to other 'mri2msh' script")
126# check if 'mri' field exists
127if not hasattr(subject, 'mri'):
128 logprint(f'\'MRI\' attribute missing in {fn_subject}. Quitting.')
129 sys.exit(1)
131mesh_folder = mesh_dict["mesh_folder"]
132if not os.path.exists(mesh_folder):
133 os.makedirs(mesh_folder)
135t1_path = subject.mri[mri_idx]['fn_mri_T1']
136t2_path = subject.mri[mri_idx]['fn_mri_T2']
137bvec = subject.mri[mri_idx]['fn_mri_DTI_bvec']
138bval = subject.mri[mri_idx]['fn_mri_DTI_bval']
139dti_path = subject.mri[mri_idx]['fn_mri_DTI']
140try:
141 dti_tensor_path = subject.mri[mri_idx]['dti_tensor_fn'] # already prepared tensor file
142except KeyError:
143 dti_tensor_path = None
145try:
146 fsdir = subject.mri[mri_idx]['freesurfer_dir']
147 if fsdir.endswith('/') or fsdir.endswith('\\'):
148 fsdir = fsdir[:-1]
149except (KeyError, AttributeError):
150 fsdir = None
152try:
153 smooth_skin = mesh_dict['smooth_skin']
154except KeyError:
155 smooth_skin = False
156if smooth_skin is None or smooth_skin == 0:
157 smooth_skin = False
158assert smooth_skin is False or 0 < smooth_skin <= 1, f"['smooth_skin'] has to be within [0,1] or None."
160if 'fn_mri_DTI_rev' in subject.mri[mri_idx].keys():
161 dti_rev_path = subject.mri[mri_idx]['fn_mri_DTI_rev']
162else:
163 dti_rev_path = ''
164if dti_rev_path is None:
165 dti_rev_path = ''
167if t2_path is None:
168 t2_path = ''
170# check if FreeSurfer env is available
171try:
172 fsaverage_dir = os.path.join(f"{os.environ['SUBJECTS_DIR']}", "fsaverage")
173except KeyError:
174 warnings.warn("Freesurfer not found, some meshing functions won't work.")
175 time.sleep(2)
177if use_fs and fsdir is None:
178 raise ValueError("use_fs is set to True but no FreeSurfer [fsdir] directory is provided.")
180# create mesh from mri if not already done
181#############################################################################################################
182if use_fs:
183 # FreeSurfer needs a subject ID so we'll use the subdirectory as the subject ID
184 fsdir_tmp, sub_id_tmp = os.path.split(fsdir)
185 if not os.path.exists(fsdir) or len(os.listdir(fsdir)) == 0: # if directory does not exist or is empty
186 fs_cmd = (f"recon-all "
187 f"-s {sub_id_tmp} "
188 f"-sd {fsdir_tmp} "
189 f"-i {subject.mri[mri_idx]['fn_mri_T1']} ")
190 if subject.mri[mri_idx]['fn_mri_T2'] is not None:
191 fs_cmd += f"-T2 {subject.mri[mri_idx]['fn_mri_T2']} "
192 fs_cmd += " -all"
193 logprint(f"Calling {fs_cmd}")
195 # Execute the recon-all command
196 ret = os.system(fs_cmd)
197 if ret == 0:
198 logprint(f"recon-all completed successfully for {subject.id}.")
199 else:
200 logprint(f"recon-all failed for {subject.id} with code {ret}.")
201 quit(ret)
202 else:
203 logprint(f"Found FreeSurfer directory {fsdir}. Skipping FreeSurfer reconstruction.")
205if approach != "charm":
206 logprint(f"Choose approach='charm' in {fn_subject} or use different script for meshing")
207 quit()
209mesh_fn = mesh_dict["fn_mesh_msh"]
210mesh_fn_hdf5 = mesh_dict['fn_mesh_hdf5']
212mesh_fn_charm = os.path.join(mesh_folder, 'm2m_' + subject.id, subject.id + ".msh")
214mesh_fn_org = mesh_fn.replace('.msh', '_org.msh')
215mesh_fn_org_hdf5 = mesh_fn_org.replace('.msh', '.hdf5')
217mesh_fn_refined = mesh_fn.replace(".msh", "_refined.msh")
218mesh_fn_refined_hdf5 = mesh_fn_refined.replace(".msh", ".hdf5")
220mesh_fn_smoothed = mesh_fn.replace('.msh', '_smoothed.msh')
221mesh_fn_smoothed_hdf5 = mesh_fn_smoothed.replace('.msh', '.hdf5')
223smooth_settings_fn = f"{mesh_folder}/m2m_{subject.id}/smooth_settings.ini"
224refinement_nii_fn = f"{mesh_folder}/m2m_{subject.id}/refinement.nii.gz"
225refinement_upsampled_nii_fn = f"{mesh_folder}/m2m_{subject.id}/refinement_upsampled.nii.gz"
227mesh_fn_unfixed = mesh_fn.replace('.msh', '_unfixed.msh')
229if not os.path.isfile(mesh_fn):
230 # open logfile
231 logprint(f"Started mri2msh_charm.py on {datetime.datetime.now().strftime('%Y-%m-%d %H:%M')}")
232 logprint("================================================")
233 logprint(f"Subject object file: {fn_subject}")
234 logprint(f"mri_idx: {mri_idx}")
235 logprint(f"mesh_idx: {mesh_idx}")
236 logprint(f"approach: {approach}")
237 logprint(f"use_fs: {use_fs}")
238 logprint(f"fs_dir: {fsdir}")
239 logprint(f"fn_mri_T1: {subject.mri[mri_idx]['fn_mri_T1']}")
240 logprint(f"fn_mri_T2: {subject.mri[mri_idx]['fn_mri_T2']}")
241 logprint(f"fn_mri_DTI_bvec: {subject.mri[mri_idx]['fn_mri_DTI_bvec']}")
242 logprint(f"fn_mri_DTI_bval: {subject.mri[mri_idx]['fn_mri_DTI_bval']}")
243 logprint(f"fn_mri_DTI: {subject.mri[mri_idx]['fn_mri_DTI']}\n")
245 # meshing -----------------------------------------------------------
246 os.chdir(mesh_folder)
247 logprint(f"Target folder: {mesh_folder}")
249 if not mesh_only:
250 # set qform to sform of T1 ------------------------------------------
251 # some scanners may yield different q- and s-form, enforce s-form here for consistency
252 img = nib.load(subject.mri[mri_idx]['fn_mri_T1'])
253 affine, code = img.get_sform(coded=True)
254 img.set_qform(affine=affine, code=code)
255 fn_t1_ext = pynibs.expio.splitext_niigz(subject.mri[mri_idx]['fn_mri_T1'])
256 subject.mri[mri_idx]['fn_mri_T1'] = fn_t1_ext[0] + "_qform" + fn_t1_ext[1]
257 nib.save(img, subject.mri[mri_idx]['fn_mri_T1'])
259 # coregister T2 to T1 -----------------------------------------------
260 # additional coregistration step improves headmodels
261 if len(t2_path) != 0:
262 logprint(f"Coregistering T2 to T1 ...")
263 cmd_charm = f"{charm_cmdline} {subject.id} {t1_path} {t2_path} " \
264 f"--registerT2 --forceqform --forcerun {debug_flag}"
265 logprint(f"Calling {cmd_charm}")
266 ret = os.system(cmd_charm)
267 if ret == 0:
268 logprint(f"Coregistering finished for {subject.id}.")
269 else:
270 logprint(f"Coregistering finished failed for {subject.id} with code {ret}")
271 quit(ret)
272 t1_path = f"{mesh_folder}/T1.nii.gz"
273 t2_path = f"{mesh_folder}/T2_reg.nii.gz"
274 os.rename(f"{mesh_folder}/m2m_{subject.id}/T1.nii.gz", t1_path)
275 os.rename(f"{mesh_folder}/m2m_{subject.id}/T2_reg.nii.gz", t2_path)
276 shutil.rmtree(f"{mesh_folder}/m2m_{subject.id}")
278 # and now set the qform again
279 img = nib.load(t1_path)
280 img.header['qform_code'] = 1
281 nib.save(img, t1_path)
283 img = nib.load(t2_path)
284 img.header['qform_code'] = 1
285 nib.save(img, t2_path)
287 # meshing -----------------------------------------------------------
288 cmd_charm = f"{charm_cmdline} {subject.id} {t1_path} {t2_path} {charm_ini_fn} {charm_cmd_params} {debug_flag}"
289 else:
290 cmd_charm = f"{charm_cmdline} {subject.id} --mesh {debug_flag}"
292 if use_fs:
293 cmd_charm += f" --fs-dir {fsdir} "
295 logprint(f"Meshing with cmd: {cmd_charm}")
296 ret = os.system(cmd_charm)
297 if ret == 0:
298 logprint(f"Meshing finished for {subject.id}.")
299 else:
300 logprint(f"Meshing finished failed for {subject.id} with code {ret}")
301 quit(ret)
303 # check mesh for degenerated elements
304 zero_tris, zero_tets, neg_tets = pynibs.util.simnibs_io.check_mesh(mesh_fn)
305 if zero_tris.size or zero_tets.size or neg_tets.size:
306 warnings.warn(f"Degenerated mesh {mesh_fn}: zero_tris: "
307 f"{zero_tris.size} zero tris, {zero_tets.size} zero_tets, {neg_tets.size} neg_tets left over.")
308 # rename .msh file to user specific filename
309 if mesh_fn_charm != mesh_fn:
310 os.rename(mesh_fn_charm, mesh_fn)
312 # create hdf5 for mesh
313 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
314 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
315 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
317 print_mesh_info(mesh_fn_hdf5)
319 # remove smoothed and refined meshes bc we created a new mesh.
320 for fn in [mesh_fn_smoothed, mesh_fn_refined, mesh_fn_smoothed_hdf5, mesh_fn_refined_hdf5]:
321 try:
322 os.unlink(fn)
323 logprint(f"Found {fn}. Removing.")
324 except FileNotFoundError:
325 pass
327elif os.path.isfile(mesh_fn):
328 logprint(f"{mesh_fn} exists. Skipping Meshing.")
329 # if not os.path.exists(mesh_fn_org):
330 # logprint(f"Orginal mesh copy not found. Creating {mesh_fn_org}.")
331 #
332 # shutil.copy(mesh_fn, mesh_fn_org)
333 # logprint(f" > Writing file: {mesh_fn_hdf5}\n")
334 # pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_org, mesh_idx=mesh_idx, skip_roi=True)
335 # pynibs.write_xdmf(mesh_fn_org_hdf5, overwrite_xdmf=True, verbose=True)
337 if not os.path.exists(mesh_fn_hdf5):
338 logprint(f"{mesh_fn_hdf5} not found.")
340 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
341 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
342 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
344if not os.path.exists(mesh_fn_org):
345 # keep a copy of the original mesh file for housekeeping reasons
346 shutil.copy(mesh_fn, mesh_fn_org)
347 logprint(f" > Writing file: {mesh_fn_org_hdf5}\n")
348 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_org, mesh_idx=mesh_idx, skip_roi=True)
349 pynibs.write_xdmf(mesh_fn_org_hdf5, overwrite_xdmf=True, verbose=True)
351# Create ROI
352############################################################################################################
353if mesh_idx in subject.roi.keys():
354 for roi_name, roi in subject.roi[mesh_idx].items():
355 time_start_roi = time.time()
357 # only run ROI things if fn_mask does not exist yet
358 if roi['fn_mask'] is not None and not os.path.exists(os.path.join(mesh_folder, roi['fn_mask'])):
360 # Default to freesurfer ROI creation to keep backwards compat
361 if roi['template'] in [None, 'fsaverage']:
362 roi_dir = os.path.join(mesh_folder, "roi", str(roi_name))
364 os.makedirs(roi_dir, exist_ok=True)
366 logprint(f"Creating ROI '{roi_name}'")
367 logprint(f"===============")
369 # write bash script with the FreeSurfer commands
370 f = open(os.path.join(roi_dir, "get_mask.sh"), "w")
371 f.write(f"export SUBJECTS_DIR='{mesh_folder}'\n")
372 f.write(f"cd {mesh_folder}\n")
373 f.write(f"rm -f fsaverage\n")
374 f.write(f"rm -f m2m_{subject.id}/surf\n")
376 f.write(f"ln -s {fsaverage_dir} 'fsaverage'\n")
377 f.write(f"ln -s 'surfaces' 'm2m_{subject.id}/surf'\n")
379 f.write(f"mri_surf2surf "
380 f"--srcsubject 'fsaverage' "
381 f"--srcsurfval '{subject.roi[mesh_idx][roi_name]['fn_mask_avg']}' "
382 f"--trgsurfval {os.path.join(mesh_folder, roi['fn_mask'])} "
383 f"--hemi {roi['hemisphere']} "
384 f"--trgsubject m2m_{subject.id} "
385 f"--trgsurfreg sphere.reg.gii "
386 f"--srcsurfreg sphere.reg\n")
387 f.close()
389 # run FreeSurfer script
390 command = f"sh {os.path.join(roi_dir, 'get_mask.sh')}"
391 logprint(f"> Running '{command}'")
392 os.popen(command).read()
393else:
394 print("No ROIs found in subject object.")
396# Refine mesh
397#####################################################
398refined = False
399if 'refinement_roi' in mesh_dict and mesh_dict['refinement_roi'] not in [None, '']:
400 refined = True
401 # create sphere that defines the refinement region
402 final_tissues_nii = f"{mesh_folder}/m2m_{subject.id}/final_tissues.nii.gz"
403 final_tissues_upsampled_nii = f"{mesh_folder}/m2m_{subject.id}/label_prep/tissue_labeling_upsampled.nii.gz"
404 refinement_settings_fn = f"{mesh_folder}/m2m_{subject.id}/refinement_size.ini"
406 if os.path.exists(mesh_fn_refined):
407 logprint(f"Refined mesh {mesh_fn_refined} exists. Skipping refinement.")
408 refine_now = False
410 elif os.path.exists(refinement_upsampled_nii_fn):
411 logprint(f"Using {refinement_upsampled_nii_fn} to refine.")
412 refine_now = True
414 else:
415 logprint("Creating new refinement (upsampled) nii.")
416 refine_now = True
417 tissue_type_2_refine = [1, 2, 3]
419 # get center of refinement sphere from roi
420 refinement_roi = mesh_dict['refinement_roi']
421 fn_roi = f"{mesh_folder}/{subject.roi[mesh_idx][refinement_roi]['fn_mask']}"
422 refined_size = mesh_dict['refinemement_element_size']
424 hem_roi_ref = subject.roi[mesh_idx][refinement_roi]['hemisphere']
425 midlayer_fn = mesh_folder + '/' + mesh_dict[f'fn_{hem_roi_ref}_midlayer']
426 roi_idx = np.where(nib.load(fn_roi).get_fdata()[:, 0, 0] != 0)[0]
427 roi_points = nib.load(midlayer_fn).agg_data()[0]
428 roi_center = np.mean(roi_points[roi_idx], axis=0)
430 with open(refinement_settings_fn, 'w') as f:
431 f.write(f"{refined_size}")
433 # find maximum distance to center
434 roi_radius = np.max(np.linalg.norm((roi_points[roi_idx] - roi_center), axis=1))
436 # create spherical refinement .nii (T1 resolution)
437 pynibs.create_refine_spherical_roi(center=roi_center,
438 radius=roi_radius,
439 target_size=refined_size,
440 final_tissues_nii=final_tissues_nii,
441 out_fn=refinement_nii_fn,
442 outside_size=0,
443 outside_factor=3,
444 verbose=True)
446 # create spherical refinement .nii (upsampled resolution)
447 pynibs.create_refine_spherical_roi(center=roi_center,
448 radius=roi_radius,
449 target_size=refined_size,
450 final_tissues_nii=final_tissues_upsampled_nii,
451 out_fn=refinement_upsampled_nii_fn,
452 outside_size=0,
453 outside_factor=3,
454 verbose=True)
456 if refine_now:
458 cmd = f'{meshmesh_cmdline} {final_tissues_upsampled_nii} ' \
459 f'{mesh_fn_refined} --sizing_field={refinement_upsampled_nii_fn} {charm_ini_fn} ' + debug_flag
461 logprint(f"Refinement: running {cmd}.")
462 os.popen(cmd).read()
464 # check mesh for degenerated elements
465 zero_tris, zero_tets, neg_tets = pynibs.util.simnibs_io.check_mesh(mesh_fn_refined)
466 if zero_tris.size or zero_tets.size or neg_tets.size:
467 warnings.warn(f"Degenerated mesh {mesh_fn_refined}:\n "
468 f"{zero_tris.size} zero tris, "
469 f"{zero_tets.size} zero_tets, "
470 f"{neg_tets.size} neg_tets identified.")
472 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
473 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_refined, mesh_idx=mesh_idx, skip_roi=True)
474 pynibs.write_xdmf(mesh_fn_refined_hdf5, overwrite_xdmf=True, verbose=True)
476 print_mesh_info(mesh_fn_refined_hdf5)
478 for fn in [mesh_fn_smoothed, mesh_fn_smoothed_hdf5, mesh_fn, mesh_fn_hdf5]:
479 try:
480 os.unlink(fn)
481 logprint(f"Found {fn}. Removing.")
482 except FileNotFoundError:
483 pass
485 shutil.copy(mesh_fn_refined, mesh_fn)
486 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
487 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
488 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
490# Mesh smoothing
491#####################################################
492if smooth_skin:
493 if os.path.exists(mesh_fn_smoothed):
494 f"f{mesh_fn_smoothed} found. Skipping smoothing."
495 else:
496 logprint(f"Smoothing skin surface with taubin filtering ({smooth_skin}).")
497 # write smooth_setting.ini file to store smooth settings
498 with open(smooth_settings_fn, 'w') as f:
499 f.write(f"{smooth_skin}")
501 smooth_mesh(mesh_fn,
502 mesh_fn_smoothed,
503 smooth=smooth_skin)
505 # check mesh for degenerated elements
506 zero_tris, zero_tets, neg_tets = pynibs.util.simnibs_io.check_mesh(mesh_fn_smoothed)
507 if zero_tris.size or zero_tets.size or neg_tets.size:
508 warnings.warn(f"Degenerated mesh {mesh_fn_smoothed}:\n "
509 f"{zero_tris.size} zero tris, "
510 f"{zero_tets.size} zero_tets, "
511 f"{neg_tets.size} neg_tets identified.")
513 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
514 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_smoothed, mesh_idx=mesh_idx, skip_roi=True)
516 pynibs.write_xdmf(mesh_fn_smoothed_hdf5, overwrite_xdmf=True, verbose=True)
518 for fn in [mesh_fn, mesh_fn_hdf5]:
519 try:
520 os.unlink(fn)
521 logprint(f"Found {fn}. Removing.")
522 except FileNotFoundError:
523 pass
525 print_mesh_info(mesh_fn_smoothed_hdf5)
527 shutil.copy(mesh_fn_smoothed, mesh_fn)
528 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
529 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
530 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
532os.chdir(mesh_folder)
534# fix mesh if neccessary
535# =================================================================================================================
536zero_tris, zero_tets, neg_tets = pynibs.util.simnibs_io.check_mesh(mesh_fn, verbose=True)
537if zero_tris.size or zero_tets.size or neg_tets.size:
538 logprint(f"Fixing degenerated elements.")
539 for fn in [mesh_fn_hdf5, mesh_fn_unfixed]:
540 try:
541 os.unlink(fn)
542 except FileNotFoundError:
543 pass
544 shutil.move(mesh_fn, mesh_fn_unfixed)
545 mesh = pynibs.util.simnibs_io.fix_mesh(mesh_fn_unfixed)
546 mesh.write(mesh_fn)
548 logprint(
549 f"Fixed {zero_tris.size} zero area triangles, {zero_tets.size} zero volume tets, {neg_tets.size}"
550 f" negative volume tets")
552 # check if new mesh is fine
553 zero_tris, zero_tets, neg_tets = pynibs.util.simnibs_io.check_mesh(mesh_fn)
554 if zero_tris.size or zero_tets.size or neg_tets.size:
555 warnings.warn(f"Couldn't fix all elements for {mesh_fn}:\n "
556 f"{zero_tris.size} zero tris, "
557 f"{zero_tets.size} zero_tets, "
558 f"{neg_tets.size} neg_tets left over.")
560# create *.hdf5 mesh file (incl. ROIs) from *.msh file
561#############################################################################################################
562logprint("Transforming mesh from .msh to .hdf5")
563logprint(f"{'=' * 36}")
565skip_msh2hdf5 = False
567if os.path.exists(mesh_fn_hdf5):
568 skip_msh2hdf5 = True
569 logprint(f" > {mesh_fn_hdf5} already exists. Checking existence of ROIs ...")
571 with h5py.File(mesh_fn_hdf5, "r") as f:
572 if mesh_idx in subject.roi.keys():
573 for roi_name in subject.roi[mesh_idx].keys():
574 try:
575 _ = f[f"roi_surface/{roi_name}"]
576 logprint(f" > roi_surface/{roi_name} exists")
578 except KeyError:
579 logprint(f" > roi_surface/{roi_name} does not exist")
580 skip_msh2hdf5 = False
581 break
582else:
583 skip_msh2hdf5 = False
585if not skip_msh2hdf5:
586 logprint(f" > Transforming .msh to .hdf5 and adding ROIs")
587 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_layer=np.logical_not(layer))
588 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
590 # create .xdmf file for visualization in paraview
591 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
592 logprint(f" > Created .xdmf file: {os.path.splitext(mesh_fn_hdf5)[0] + '.xdmf'}")
594else:
595 logprint(" > Skipping transformation from .msh to .hdf5 (already exists and all ROIs are included)")
597roi = pynibs.load_roi_surface_obj_from_hdf5(fname=subject.mesh[mesh_idx]['fn_mesh_hdf5'])
598if roi is not None:
599 for roi_name in subject.roi[mesh_idx].keys():
600 fn_geo_roi = os.path.join(mesh_folder, "roi", roi_name, "geo.hdf5")
602 if not os.path.exists(fn_geo_roi):
603 logprint(f" > Creating geo.hdf5 for ROI #{roi_name}: {fn_geo_roi}")
605 if not os.path.exists(os.path.split(fn_geo_roi)[0]):
606 os.makedirs(os.path.split(fn_geo_roi)[0])
608 # midlayer .geo file
609 pynibs.write_geo_hdf5_surf(out_fn=fn_geo_roi,
610 points=roi[roi_name].node_coord_mid,
611 con=roi[roi_name].node_number_list,
612 replace=True,
613 hdf5_path='/mesh')
614 pynibs.write_xdmf(fn_geo_roi)
616# determine anisotropic conductivity tensors
617#############################################################################################################
618d2c_folder = os.path.join(mesh_folder, f'm2m_{subject.id}', 'dMRI_prep')
619create_conductivity_tensors = bvec is not None and \
620 bval is not None and \
621 dti_path is not None \
622 and not (os.path.exists(os.path.join(mesh_dict["mesh_folder"],
623 mesh_dict['fn_tensor_vn']))
624 and os.path.exists(os.path.join(d2c_folder,
625 "first_ev_for_check.hdf5"))
626 )
627register_existing_tensor = dti_tensor_path is not None \
628 and not (os.path.exists(os.path.join(mesh_dict["mesh_folder"],
629 mesh_dict['fn_tensor_vn']))
630 and os.path.exists(os.path.join(d2c_folder,
631 "first_ev_for_check.hdf5"))
632 )
634if create_conductivity_tensors and not register_existing_tensor:
635 time_start_dwi = time.time()
636 tensor_fn = os.path.join(mesh_folder, f'm2m_{subject.id}', "dMRI_prep", "first_ev_for_check.msh")
637 tensor_fn_hdf5 = tensor_fn.replace('.msh', 'hdf5')
638 if not os.path.exists(os.path.join(mesh_folder, mesh_dict["fn_tensor_vn"])):
639 try:
640 os.remove(tensor_fn_hdf5)
641 except FileNotFoundError:
642 pass
644 # fix .bvals to 0:
645 assert os.path.exists(subject.mri[mri_idx]['fn_mri_DTI_bval']), \
646 f"Can't find {subject.mri[mri_idx]['fn_mri_DTI_bval']}"
648 with open(subject.mri[mri_idx]['fn_mri_DTI_bval'], 'r+') as f:
649 f.seek(0)
650 bvals = f.read()
652 bvals = bvals.replace('\n', '')
654 if bvals[-1] == " ":
655 bvals = bvals[:-1]
657 bvals = np.array(bvals.split(" ")).astype(int)
659 if 0 not in bvals:
660 min_val = np.min(bvals)
661 assert min_val <= 10
662 logprint(f"No b0 found in {subject.mri[mri_idx]['fn_mri_DTI_bval']}. Changing {min_val} to 0.")
663 bvals[bvals == min_val] = 0
664 bvals = " ".join(bvals.astype(str).tolist())
665 shutil.copy(subject.mri[mri_idx]['fn_mri_DTI_bval'], subject.mri[mri_idx]['fn_mri_DTI_bval'] + '.org')
666 f.seek(0)
667 f.write(bvals)
669 logprint("Charm version of dwi2cond")
671 command = f"{dwi2cond_cmdline} " \
672 f"--all " \
673 f"--eddy " \
674 f"--readout={subject.mri[mri_idx]['dti_readout_time']}" \
675 f" --phasedir={subject.mri[mri_idx]['dti_phase_direction']} " \
676 f"{subject.id} {dti_path} {bval} {bvec} {dti_rev_path}"
678 logprint("Creating anisotropic conductivity tensors from DWI data")
679 logprint("=======================================================")
680 logprint(f"Running {command}")
682 os.putenv('SUBJECTS_DIR', os.getcwd())
684 os.popen(command).read()
685 logprint(f" > Created files in folder {os.path.join(mesh_folder, f'm2m_{subject.id}', 'dMRI_prep')}\n")
687 elif register_existing_tensor:
688 logprint(f"Using preprocess tensor file {dti_tensor_path} to create conductivity tensors.")
689 time_start_dwi = time.time()
690 tensor_fn = os.path.join(mesh_folder, f'm2m_{subject.id}', "dMRI_prep", "first_ev_for_check.msh")
691 tensor_fn_hdf5 = tensor_fn.replace('.msh', 'hdf5')
692 if not os.path.exists(os.path.join(mesh_folder, mesh_dict["fn_tensor_vn"])):
693 try:
694 os.remove(tensor_fn_hdf5)
695 except FileNotFoundError:
696 pass
698 command = f"{dwi2cond_cmdline} " \
699 f"--all " \
700 f"{subject.id} {dti_tensor_path} "
702 logprint("Creating anisotropic conductivity tensors from DWI data")
703 logprint("=======================================================")
704 logprint(f"Running {command}")
706 # is this also necessary here? Not sure
707 os.putenv('SUBJECTS_DIR', os.getcwd())
709 os.popen(command).read()
710 logprint(f" > Created files in folder {os.path.join(mesh_folder, f'm2m_{subject.id}', 'dMRI_prep')}\n")
711 else:
712 logprint(f" > Anisotropy files already exist ... skipping to create "
713 f"{mesh_dict['fn_tensor_vn']}")
715 time_stop_dwi = time.time()
716 time_dwi = time_stop_dwi - time_start_dwi
718 if not os.path.exists(tensor_fn_hdf5):
719 # write .hdf5 tensor file
720 pynibs.msh2hdf5(subject=subject,
721 mesh_idx=mesh_idx,
722 fn_msh=tensor_fn,
723 skip_roi=True,
724 include_data=True)
726logprint(f"Time: Creating mesh: {time_mesh:6.2f} secs")
727if create_conductivity_tensors:
728 logprint(f"Time: Creating FA: {time_dwi:6.2f} secs")
729logprint("============ DONE ============")