Coverage for /builds/tms-localization/papers/tmsloc_proto/scripts/02_make_msh_from_mri_simnibs4.py: 69%
392 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-17 07:31 +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 pynibs.util.simnibs import smooth_mesh
16# Get command line arguments
17parser = argparse.ArgumentParser(description='Creates a FEM mesh from MRI data')
18parser.add_argument('-s', '--fn_subject', help='Path to *.hdf5-file of subject', required=True)
19parser.add_argument('--charm_ini', help='charm.ini file to use.',
20 default=f"{os.path.dirname(os.path.realpath(__file__))}/configs/charm_tms410.ini")
21parser.add_argument('-m', '--mesh_idx', help='The mesh name from subject object', required=True)
22parser.add_argument('-l', '--layer', help='Create cortical layer', action='store_true', required=False, default=False)
23parser.add_argument('--charm_cmd_params', help='Additional parameters to charm, e.g. forcerun or mesh', required=False, default='')
24parser.add_argument('--mesh', help='Skip segmentation and only mesh headmodel', action='store_true', required=False, default=False)
25parser.add_argument('--debug', help='Run charm in debug mode saving additional intermediate output', action='store_true', required=False, default=False)
27args = parser.parse_args()
28charm_cmdline = "charm"
29dwi2cond_cmdline = "dwi2cond"
30if os.system(f"which {dwi2cond_cmdline}") == 256:
31 simnibs_dir = os.path.split(simnibs.__file__)[0]
32 dwi2cond_cmdline = f"{simnibs_dir}/external/dwi2cond"
33 assert os.system(f"which {dwi2cond_cmdline}") == 0, f"Cannot find dwi2cond executable."
34meshmesh_cmdline = "meshmesh"
35layer = args.layer
36mesh_only = args.mesh
38debug_flag = ""
39if args.debug:
40 debug_flag = "--debug"
43def logprint(text):
44 """Just a wrapper to print and write text to a logfile."""
45 with open(log_file, "a") as logfile:
46 logfile.write(f'{text}\n')
47 print(text)
49fn_subject = os.path.abspath(args.fn_subject)
50subject = pynibs.load_subject(fn_subject)
51mesh_idx = args.mesh_idx
52charm_ini_fn = args.charm_ini
54assert mesh_idx in subject.mesh, f"Mesh '{mesh_idx}' not found in {args.fn_subject}"
55mesh_folder = subject.mesh[mesh_idx]["mesh_folder"]
56log_file = os.path.join(mesh_folder, f"info_{datetime.datetime.now():%Y%m%d_%H%M%S}.log")
57charm_cmd_params = args.charm_cmd_params
58if charm_cmd_params != '' and not charm_cmd_params.startswith('--'):
59 charm_cmd_params = '--' + charm_cmd_params
62def print_mesh_info(fn):
63 logprint("- " * 32)
64 logprint(f"Mesh information: {fn}")
65 pref = os.path.basename(fn)
66 with h5py.File(fn, 'r') as f: # f = h5py.File(mesh_fn_smoothed_hdf5,'r')
67 logprint(f"{pref}: total nodes: {f['/mesh/nodes/node_number'].shape[0]: >7}")
68 logprint(f"{pref}: total elms: {f['/mesh/elm/elm_number'].shape[0]: >7}")
69 logprint(f"{pref}: total tris: {np.sum(f['/mesh/elm/elm_type'][:] == 2): >7}")
70 logprint(f"{pref}: total tets: {np.sum(f['/mesh/elm/elm_type'][:] == 4): >7}")
71 logprint("- " * 32)
72 for t in np.unique(f['/mesh/elm/tet_tissue_type']):
73 logprint(f"{pref}: tet tissue type {t:0>4}: {np.sum(f['/mesh/elm/tet_tissue_type'] == t): >7}")
74 logprint("- " * 32)
75 for t in np.unique(f['/mesh/elm/tri_tissue_type']):
76 logprint(f"{pref}: tri tissue type {t:0>4}: {np.sum(f['/mesh/elm/tri_tissue_type'] == t): >7}")
77 logprint("- " * 32)
80time_mesh = 0
81time_roi = 0
82time_dwi = 0
84# logprint simulation info
85# ================================================================================
86os.makedirs(os.path.split(log_file)[0], exist_ok=True)
87logprint("-" * 64)
88logprint(f"{parser.description}")
89args_dict = vars(args)
90for key in args_dict.keys():
91 logprint(f"{key: >15}: {args_dict[key]}")
92logprint("-" * 64)
94if charm_ini_fn:
95 # this will always be appended to all charm, meshmesh calls
96 charm_ini_fn = f" --usesettings={os.path.abspath(charm_ini_fn)}"
97else:
98 charm_ini_fn = ''
100# load subject info
101mesh_dict = subject.mesh[mesh_idx]
102mri_idx = mesh_dict["mri_idx"]
104# set approach specific options
105approach = mesh_dict["approach"]
107if approach == "charm":
108 vertex_density = None
109 numvertices = None
110else:
111 raise NotImplementedError(f"Approach {approach} not implemented. Use 'charm' or switch to other 'mri2msh' script")
113# check if 'mri' field exists
114if not hasattr(subject, 'mri'):
115 logprint(f'\'MRI\' attribute missing in {fn_subject}. Quitting.')
116 sys.exit(1)
118mesh_folder = mesh_dict["mesh_folder"]
119if not os.path.exists(mesh_folder):
120 os.makedirs(mesh_folder)
122t1_path = subject.mri[mri_idx]['fn_mri_T1']
123t2_path = subject.mri[mri_idx]['fn_mri_T2']
124bvec = subject.mri[mri_idx]['fn_mri_DTI_bvec']
125bval = subject.mri[mri_idx]['fn_mri_DTI_bval']
126dti_path = subject.mri[mri_idx]['fn_mri_DTI']
128try:
129 fsdir = subject.mri[mri_idx]['freesurfer_dir']
130except KeyError:
131 fsdir = None
132if fsdir == '':
133 fsdir = None
135if fsdir is not None:
136 assert os.path.exists(fsdir), f"Cannot find freesurfer directory for this subject {fsdir}."
138try:
139 smooth_skin = mesh_dict['smooth_skin']
140except KeyError:
141 smooth_skin = False
142if smooth_skin is None or smooth_skin == 0:
143 smooth_skin = False
144assert smooth_skin == False or 0 < smooth_skin <= 1, f"['smooth_skin'] has to be within [0,1] or None."
146if 'fn_mri_DTI_rev' in subject.mri[mri_idx].keys():
147 dti_rev_path = subject.mri[mri_idx]['fn_mri_DTI_rev']
148else:
149 dti_rev_path = ''
150if dti_rev_path is None:
151 dti_rev_path = ''
153if t2_path is None:
154 t2_path = ''
156# check if FreeSurfer env is available
157try:
158 fsaverage_dir = os.path.join(f"{os.environ['SUBJECTS_DIR']}", "fsaverage")
159except KeyError:
160 warnings.warn("Freesurfer not found.")
162# create mesh from mri if not already done
163#############################################################################################################
164if approach != "charm":
165 logprint(f"Choose approach='charm' in {fn_subject} or use different script for meshing")
166 quit()
168mesh_fn = mesh_dict["fn_mesh_msh"]
169mesh_fn_hdf5 = mesh_dict['fn_mesh_hdf5']
171mesh_fn_charm = os.path.join(mesh_folder, 'm2m_' + subject.id, subject.id + ".msh")
173mesh_fn_org = mesh_fn.replace('.msh', '_org.msh')
174mesh_fn_org_hdf5 = mesh_fn_org.replace('.msh', '.hdf5')
176mesh_fn_refined = mesh_fn.replace(".msh", "_refined.msh")
177mesh_fn_refined_hdf5 = mesh_fn_refined.replace(".msh", ".hdf5")
179mesh_fn_smoothed = mesh_fn.replace('.msh', '_smoothed.msh')
180mesh_fn_smoothed_hdf5 = mesh_fn_smoothed.replace('.msh', '.hdf5')
182smooth_settings_fn = f"{mesh_folder}/m2m_{subject.id}/smooth_settings.ini"
183refinement_nii_fn = f"{mesh_folder}/m2m_{subject.id}/refinement.nii.gz"
184refinement_upsampled_nii_fn = f"{mesh_folder}/m2m_{subject.id}/refinement_upsampled.nii.gz"
186mesh_fn_unfixed = mesh_fn.replace('.msh', '_unfixed.msh')
188if not os.path.isfile(mesh_fn):
189 # open logfile
190 logprint(f"Started mri2msh_charm.py on {datetime.datetime.now().strftime('%Y-%m-%d %H:%M')}")
191 logprint("================================================")
192 logprint(f"Subject object file: {fn_subject}")
193 logprint(f"mri_idx: {mri_idx}")
194 logprint(f"mesh_idx: {mesh_idx}")
195 logprint(f"approach: {approach}")
196 logprint(f"fn_mri_T1: {subject.mri[mri_idx]['fn_mri_T1']}")
197 logprint(f"fn_mri_T2: {subject.mri[mri_idx]['fn_mri_T2']}")
198 logprint(f"fn_mri_DTI_bvec: {subject.mri[mri_idx]['fn_mri_DTI_bvec']}")
199 logprint(f"fn_mri_DTI_bval: {subject.mri[mri_idx]['fn_mri_DTI_bval']}")
200 logprint(f"fn_mri_DTI: {subject.mri[mri_idx]['fn_mri_DTI']}\n")
202 # meshing -----------------------------------------------------------
203 os.chdir(mesh_folder)
204 logprint(f"Target folder: {mesh_folder}")
206 if not mesh_only:
207 # set qform to sform of T1 ------------------------------------------
208 # some scanners may yield different q- and s-form, enforce s-form here for consistency
209 img = nib.load(subject.mri[mri_idx]['fn_mri_T1'])
210 affine, code = img.get_sform(coded=True)
211 img.set_qform(affine=affine, code=code)
212 fn_t1_ext = pynibs.splitext_niigz(subject.mri[mri_idx]['fn_mri_T1'])
213 subject.mri[mri_idx]['fn_mri_T1'] = fn_t1_ext[0] + "_qform" + fn_t1_ext[1]
214 nib.save(img, subject.mri[mri_idx]['fn_mri_T1'])
216 # coregister T2 to T1 -----------------------------------------------
217 # additional coregistration step improves headmodels
218 if len(t2_path) != 0:
219 logprint(f"Coregistering T2 to T1 ...")
220 cmd_charm = f"{charm_cmdline} {subject.id} {t1_path} {t2_path} " \
221 f"--registerT2 --forceqform --forcerun {debug_flag}"
222 os.system(cmd_charm)
223 logprint(f"Coregistering finished.")
224 t1_path = f"{mesh_folder}/T1.nii.gz"
225 t2_path = f"{mesh_folder}/T2_reg.nii.gz"
226 os.rename(f"{mesh_folder}/m2m_{subject.id}/T1.nii.gz", t1_path)
227 os.rename(f"{mesh_folder}/m2m_{subject.id}/T2_reg.nii.gz", t2_path)
228 shutil.rmtree(f"{mesh_folder}/m2m_{subject.id}")
230 # meshing -----------------------------------------------------------
231 cmd_charm = f"{charm_cmdline} {subject.id} {t1_path} {t2_path} {charm_ini_fn} {charm_cmd_params} {debug_flag}"
232 else:
233 cmd_charm = f"{charm_cmdline} {subject.id} --mesh {debug_flag}"
234 from simnibs import __version__ as sim_v
235 if sim_v.startswith('4.1'):
236 cmd_charm += f" --fsdir {fsdir} " if fsdir is not None else ''
237 logprint(f"Meshing with cmd: {cmd_charm}.")
238 os.system(cmd_charm)
239 logprint(f"Meshing finished: {mesh_fn}.")
241 # check mesh for degenerated elements
242 zero_tris, zero_tets, neg_tets = pynibs.check_mesh(mesh_fn)
243 if zero_tris.size or zero_tets.size or neg_tets.size:
244 warnings.warn(f"Degenerated mesh {mesh_fn}: zero_tris: "
245 f"{zero_tris.size} zero tris, {zero_tets.size} zero_tets, {neg_tets.size} neg_tets left over.")
246 # rename .msh file to user specific filename
247 if mesh_fn_charm != mesh_fn:
248 os.rename(mesh_fn_charm, mesh_fn)
250 # create hdf5 for mesh
251 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
252 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
253 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
255 print_mesh_info(mesh_fn_hdf5)
257 # remove smoothed and refined meshes bc we created a new mesh.
258 for fn in [mesh_fn_smoothed, mesh_fn_refined, mesh_fn_smoothed_hdf5, mesh_fn_refined_hdf5]:
259 try:
260 os.unlink(fn)
261 logprint(f"Found {fn}. Removing.")
262 except FileNotFoundError:
263 pass
265elif os.path.isfile(mesh_fn):
266 logprint(f"{mesh_fn} exists. Skipping Meshing.")
267 # if not os.path.exists(mesh_fn_org):
268 # logprint(f"Orginal mesh copy not found. Creating {mesh_fn_org}.")
269 #
270 # shutil.copy(mesh_fn, mesh_fn_org)
271 # logprint(f" > Writing file: {mesh_fn_hdf5}\n")
272 # pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_org, mesh_idx=mesh_idx, skip_roi=True)
273 # pynibs.write_xdmf(mesh_fn_org_hdf5, overwrite_xdmf=True, verbose=True)
275 if not os.path.exists(mesh_fn_hdf5):
276 logprint(f"{mesh_fn_hdf5} not found.")
278 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
279 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
280 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
282if not os.path.exists(mesh_fn_org):
283 # keep a copy of the original mesh file for housekeeping reasons
284 shutil.copy(mesh_fn, mesh_fn_org)
285 logprint(f" > Writing file: {mesh_fn_org_hdf5}\n")
286 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_org, mesh_idx=mesh_idx, skip_roi=True)
287 pynibs.write_xdmf(mesh_fn_org_hdf5, overwrite_xdmf=True, verbose=True)
289# Create ROI
290############################################################################################################
291if mesh_idx in subject.roi.keys():
292 for roi_name, roi in subject.roi[mesh_idx].items():
293 time_start_roi = time.time()
295 # only run ROI things if fn_mask does not exist yet
296 if roi['fn_mask'] is not None and not os.path.exists(os.path.join(mesh_folder, roi['fn_mask'])):
298 # Default to freesurfer ROI creation to keep backwards compat
299 if roi['template'] in [None, 'fsaverage']:
300 roi_dir = os.path.join(mesh_folder, "roi", str(roi_name))
302 os.makedirs(roi_dir, exist_ok=True)
304 logprint(f"Creating ROI '{roi_name}'")
305 logprint(f"===============")
307 # write bash script with the FreeSurfer commands
308 f = open(os.path.join(roi_dir, "get_mask.sh"), "w")
309 f.write(f"export SUBJECTS_DIR='{mesh_folder}'\n")
310 f.write(f"cd {mesh_folder}\n")
311 f.write(f"rm -f fsaverage\n")
312 f.write(f"rm -f m2m_{subject.id}/surf\n")
314 f.write(f"ln -s {fsaverage_dir} 'fsaverage'\n")
315 f.write(f"ln -s 'surfaces' 'm2m_{subject.id}/surf'\n")
316 # f.write(
317 # f"mri_surf2surf --srcsubject 'fsaverage' "
318 # f"--srcsurfval '{subject.roi[mesh_idx][roi_name]['fn_mask_avg']}' "
319 # f"--trgsubject 'fs_'{subject.id} "
320 # f"--trgsurfval {os.path.join(mesh_folder, roi['fn_mask'])} "
321 # f"--hemi {roi['hemisphere']}\n")
322 f.write(f"mri_surf2surf "
323 f"--srcsubject 'fsaverage' "
324 f"--srcsurfval '{subject.roi[mesh_idx][roi_name]['fn_mask_avg']}' "
325 f"--trgsurfval {os.path.join(mesh_folder, roi['fn_mask'])} "
326 f"--hemi {roi['hemisphere']} "
327 f"--trgsubject m2m_{subject.id} "
328 f"--trgsurfreg sphere.reg.gii "
329 f"--srcsurfreg sphere.reg\n")
330 f.close()
332 # run FreeSurfer script
333 command = f"sh {os.path.join(roi_dir, 'get_mask.sh')}"
334 logprint(f"> Running '{command}'")
335 os.popen(command).read()
336else:
337 print("No ROIs found in subject object.")
339# Refine mesh
340#####################################################
341refined = False
342if 'refinement_roi' in mesh_dict and mesh_dict['refinement_roi'] not in [None, '']:
343 refined = True
344 # create sphere that defines the refinement region
345 final_tissues_nii = f"{mesh_folder}/m2m_{subject.id}/final_tissues.nii.gz"
346 final_tissues_upsampled_nii = f"{mesh_folder}/m2m_{subject.id}/label_prep/tissue_labeling_upsampled.nii.gz"
347 refinement_settings_fn = f"{mesh_folder}/m2m_{subject.id}/refinement_size.ini"
349 if os.path.exists(mesh_fn_refined):
350 logprint(f"Refined mesh {mesh_fn_refined} exists. Skipping refinement.")
351 refine_now = False
353 elif os.path.exists(refinement_upsampled_nii_fn):
354 logprint(f"Using {refinement_upsampled_nii_fn} to refine.")
355 refine_now = True
357 else:
358 logprint("Creating new refinement (upsampled) nii.")
359 refine_now = True
360 tissue_type_2_refine = [1, 2, 3]
362 # get center of refinement sphere from roi
363 refinement_roi = mesh_dict['refinement_roi']
364 fn_roi = f"{mesh_folder}/{subject.roi[mesh_idx][refinement_roi]['fn_mask']}"
365 refined_size = mesh_dict['refinemement_element_size']
367 hem_roi_ref = subject.roi[mesh_idx][refinement_roi]['hemisphere']
368 midlayer_fn = mesh_folder + '/' + mesh_dict[f'fn_{hem_roi_ref}_midlayer']
369 roi_idx = np.where(nib.load(fn_roi).get_fdata()[:, 0, 0] != 0)[0]
370 roi_points = nib.load(midlayer_fn).agg_data()[0]
371 roi_center = np.mean(roi_points[roi_idx], axis=0)
373 with open(refinement_settings_fn, 'w') as f:
374 f.write(f"{refined_size}")
376 # find maximum distance to center
377 roi_radius = np.max(np.linalg.norm((roi_points[roi_idx] - roi_center), axis=1))
379 # create spherical refinement .nii (T1 resolution)
380 pynibs.create_refine_spherical_roi(center=roi_center,
381 radius=roi_radius,
382 target_size=refined_size,
383 final_tissues_nii=final_tissues_nii,
384 out_fn=refinement_nii_fn,
385 outside_size=0,
386 outside_factor=3,
387 verbose=True)
389 # create spherical refinement .nii (upsampled resolution)
390 pynibs.create_refine_spherical_roi(center=roi_center,
391 radius=roi_radius,
392 target_size=refined_size,
393 final_tissues_nii=final_tissues_upsampled_nii,
394 out_fn=refinement_upsampled_nii_fn,
395 outside_size=0,
396 outside_factor=3,
397 verbose=True)
399 if refine_now:
401 cmd = f'{meshmesh_cmdline} {final_tissues_upsampled_nii} ' \
402 f'{mesh_fn_refined} --sizing_field={refinement_upsampled_nii_fn} {charm_ini_fn} ' + debug_flag
404 logprint(f"Refinement: running {cmd}.")
405 os.popen(cmd).read()
407 # check mesh for degenerated elements
408 zero_tris, zero_tets, neg_tets = pynibs.check_mesh(mesh_fn_refined)
409 if zero_tris.size or zero_tets.size or neg_tets.size:
410 warnings.warn(f"Degenerated mesh {mesh_fn_refined}:\n "
411 f"{zero_tris.size} zero tris, "
412 f"{zero_tets.size} zero_tets, "
413 f"{neg_tets.size} neg_tets identified.")
415 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
416 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_refined, mesh_idx=mesh_idx, skip_roi=True)
417 pynibs.write_xdmf(mesh_fn_refined_hdf5, overwrite_xdmf=True, verbose=True)
419 print_mesh_info(mesh_fn_refined_hdf5)
421 for fn in [mesh_fn_smoothed, mesh_fn_smoothed_hdf5, mesh_fn, mesh_fn_hdf5]:
422 try:
423 os.unlink(fn)
424 logprint(f"Found {fn}. Removing.")
425 except FileNotFoundError:
426 pass
428 shutil.copy(mesh_fn_refined, mesh_fn)
429 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
430 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
431 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
433# Mesh smoothing
434#####################################################
435if smooth_skin:
436 if os.path.exists(mesh_fn_smoothed):
437 f"f{mesh_fn_smoothed} found. Skipping smoothing."
438 else:
439 logprint(f"Smoothing skin surface with taubin filtering ({smooth_skin}).")
440 # write smooth_setting.ini file to store smooth settings
441 with open(smooth_settings_fn, 'w') as f:
442 f.write(f"{smooth_skin}")
444 smooth_mesh(mesh_fn,
445 mesh_fn_smoothed,
446 smooth=smooth_skin)
448 # check mesh for degenerated elements
449 zero_tris, zero_tets, neg_tets = pynibs.check_mesh(mesh_fn_smoothed)
450 if zero_tris.size or zero_tets.size or neg_tets.size:
451 warnings.warn(f"Degenerated mesh {mesh_fn_smoothed}:\n "
452 f"{zero_tris.size} zero tris, "
453 f"{zero_tets.size} zero_tets, "
454 f"{neg_tets.size} neg_tets identified.")
456 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
457 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn_smoothed, mesh_idx=mesh_idx, skip_roi=True)
459 pynibs.write_xdmf(mesh_fn_smoothed_hdf5, overwrite_xdmf=True, verbose=True)
461 for fn in [mesh_fn, mesh_fn_hdf5]:
462 try:
463 os.unlink(fn)
464 logprint(f"Found {fn}. Removing.")
465 except FileNotFoundError:
466 pass
468 print_mesh_info(mesh_fn_smoothed_hdf5)
470 shutil.copy(mesh_fn_smoothed, mesh_fn)
471 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
472 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_roi=True)
473 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
475os.chdir(mesh_folder)
477# fix mesh if neccessary
478# =================================================================================================================
479zero_tris, zero_tets, neg_tets = pynibs.check_mesh(mesh_fn, verbose=True)
480if zero_tris.size or zero_tets.size or neg_tets.size:
481 logprint(f"Fixing degenerated elements.")
482 for fn in [mesh_fn_hdf5, mesh_fn_unfixed]:
483 try:
484 os.unlink(fn)
485 except FileNotFoundError:
486 pass
487 shutil.move(mesh_fn, mesh_fn_unfixed)
488 mesh = pynibs.fix_mesh(mesh_fn_unfixed)
489 mesh.write(mesh_fn)
491 logprint(
492 f"Fixed {zero_tris.size} zero area triangles, {zero_tets.size} zero volume tets, {neg_tets.size} negative volume tets")
494 # check if new mesh is fine
495 zero_tris, zero_tets, neg_tets = pynibs.check_mesh(mesh_fn)
496 if zero_tris.size or zero_tets.size or neg_tets.size:
497 warnings.warn(f"Couldn't fix all elements for {mesh_fn}:\n "
498 f"{zero_tris.size} zero tris, "
499 f"{zero_tets.size} zero_tets, "
500 f"{neg_tets.size} neg_tets left over.")
502# create *.hdf5 mesh file (incl. ROIs) from *.msh file
503#############################################################################################################
504logprint("Transforming mesh from .msh to .hdf5")
505logprint(f"{'=' * 36}")
507skip_msh2hdf5 = False
509if os.path.exists(mesh_fn_hdf5):
510 skip_msh2hdf5 = True
511 logprint(f" > {mesh_fn_hdf5} already exists. Checking existence of ROIs ...")
513 with h5py.File(mesh_fn_hdf5, "r") as f:
514 if mesh_idx in subject.roi.keys():
515 for roi_name in subject.roi[mesh_idx].keys():
516 try:
517 _ = f[f"roi_surface/{roi_name}"]
518 logprint(f" > roi_surface/{roi_name} exists")
520 except KeyError:
521 logprint(f" > roi_surface/{roi_name} does not exist")
522 skip_msh2hdf5 = False
523 break
524else:
525 skip_msh2hdf5 = False
527if not skip_msh2hdf5:
528 logprint(f" > Transforming .msh to .hdf5 and adding ROIs")
529 pynibs.msh2hdf5(subject=subject, fn_msh=mesh_fn, mesh_idx=mesh_idx, skip_layer=np.logical_not(layer))
530 logprint(f" > Writing file: {mesh_fn_hdf5}\n")
532 # create .xdmf file for visualization in paraview
533 pynibs.write_xdmf(mesh_fn_hdf5, overwrite_xdmf=True, verbose=True)
534 logprint(f" > Created .xdmf file: {os.path.splitext(mesh_fn_hdf5)[0] + '.xdmf'}")
536else:
537 logprint(" > Skipping transformation from .msh to .hdf5 (already exists and all ROIs are included)")
539roi = pynibs.load_roi_surface_obj_from_hdf5(fname=subject.mesh[mesh_idx]['fn_mesh_hdf5'])
540if roi is not None:
541 for roi_name in subject.roi[mesh_idx].keys():
542 fn_geo_roi = os.path.join(mesh_folder, "roi", roi_name, "geo.hdf5")
544 if not os.path.exists(fn_geo_roi):
545 logprint(f" > Creating geo.hdf5 for ROI #{roi_name}: {fn_geo_roi}")
547 if not os.path.exists(os.path.split(fn_geo_roi)[0]):
548 os.makedirs(os.path.split(fn_geo_roi)[0])
550 # midlayer .geo file
551 pynibs.write_geo_hdf5_surf(out_fn=fn_geo_roi,
552 points=roi[roi_name].node_coord_mid,
553 con=roi[roi_name].node_number_list,
554 replace=True,
555 hdf5_path='/mesh')
556 pynibs.write_xdmf(fn_geo_roi)
558# determine anisotropic conductivity tensors
559#############################################################################################################
560d2c_folder = os.path.join(mesh_folder, f'm2m_{subject.id}', 'dMRI_prep')
561create_conductivity_tensors = bvec is not None and \
562 bval is not None and \
563 dti_path is not None \
564 and not (os.path.exists(os.path.join(mesh_dict["mesh_folder"],
565 mesh_dict['fn_tensor_vn']))
566 and os.path.exists(os.path.join(d2c_folder,
567 "first_ev_for_check.hdf5"))
568 )
570if create_conductivity_tensors:
571 time_start_dwi = time.time()
572 tensor_fn = os.path.join(mesh_folder, 'm2m_' + str(subject.id), "dMRI_prep", "first_ev_for_check.msh")
573 tensor_fn_hdf5 = tensor_fn.replace('.msh', 'hdf5')
574 if not os.path.exists(os.path.join(mesh_folder, mesh_dict["fn_tensor_vn"])):
575 try:
576 os.remove(tensor_fn_hdf5)
577 except FileNotFoundError:
578 pass
580 # fix .bvals to 0:
581 assert os.path.exists(subject.mri[mri_idx]['fn_mri_DTI_bval']), \
582 f"Can't find {subject.mri[mri_idx]['fn_mri_DTI_bval']}"
584 with open(subject.mri[mri_idx]['fn_mri_DTI_bval'], 'r+') as f:
585 f.seek(0)
586 bvals = f.read()
588 bvals = bvals.replace('\n', '')
590 if bvals[-1] == " ":
591 bvals = bvals[:-1]
593 bvals = np.array(bvals.split(" ")).astype(int)
595 if 0 not in bvals:
596 min_val = np.min(bvals)
597 assert min_val <= 10
598 logprint(f"No b0 found in {subject.mri[mri_idx]['fn_mri_DTI_bval']}. Changing {min_val} to 0.")
599 bvals[bvals == min_val] = 0
600 bvals = " ".join(bvals.astype(str).tolist())
601 shutil.copy(subject.mri[mri_idx]['fn_mri_DTI_bval'], subject.mri[mri_idx]['fn_mri_DTI_bval'] + '.org')
602 f.seek(0)
603 f.write(bvals)
605 logprint("Charm version of dwi2cond")
607 command = f"{dwi2cond_cmdline} " \
608 f"--all " \
609 f"--eddy " \
610 f"--readout={subject.mri[mri_idx]['dti_readout_time']}" \
611 f" --phasedir={subject.mri[mri_idx]['dti_phase_direction']} " \
612 f"{subject.id} {dti_path} {bval} {bvec} {dti_rev_path}"
614 logprint("Creating anisotropic conductivity tensors from DWI data")
615 logprint("=======================================================")
616 logprint(f"Running {command}")
618 os.putenv('SUBJECTS_DIR', os.getcwd())
620 os.popen(command).read()
621 logprint(f" > Created files in folder {os.path.join(mesh_folder, 'm2m_' + str(subject.id), 'dMRI_prep')}\n")
623 else:
624 logprint(f" > Anisotropy files already exist ... skipping to create "
625 f"{mesh_dict['fn_tensor_vn']}")
627 time_stop_dwi = time.time()
628 time_dwi = time_stop_dwi - time_start_dwi
630 if not os.path.exists(tensor_fn_hdf5):
631 # write .hdf5 tensor file
632 pynibs.msh2hdf5(subject=subject,
633 mesh_idx=mesh_idx,
634 fn_msh=tensor_fn,
635 skip_roi=True,
636 include_data=True)
638logprint(f"Time: Creating mesh: {time_mesh:6.2f} secs")
639if create_conductivity_tensors:
640 logprint(f"Time: Creating FA: {time_dwi:6.2f} secs")
641logprint("============ DONE ============")