The shape of the T1 image and the shape of the electric field are different
CAD Modeling
18
Posts
2
Posters
228
Views
2
Watching
-
You don't have the same resolution (grid spacing) in the T1 image and the efield...!
To visualize the two "fields" on the same grid, drag the T1w image to the analysis explorer, and interpolate it on the e-field grid.The interpolator has two inputs:
- T1w image (input field)
- E-field (target grid)
Here is a snippet to show this could be done in Python
import s4l_v1.analysis as analysis import s4l_v1.document as document import s4l_v1.model as model # Add a new ModelToGridFilter inputs = [] model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs) model_to_grid_filter.Name = "Image" model_to_grid_filter.Entity = model.AllEntities()["T1w-image"] # The model entity corresponding to the T1w-image model_to_grid_filter.UpdateAttributes() # Add a new SimulationExtractor simulation = document.AllSimulations["EM"] simulation_extractor = simulation.Results() # Add a new EmSensorExtractor em_sensor_extractor = simulation_extractor["Overall Field"] em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All" # Add a new FieldInterpolationFilter inputs = [model_to_grid_filter.Outputs[""], em_sensor_extractor.Outputs["EM E(x,y,z,f0)"]] field_interpolation_filter = analysis.core.FieldInterpolationFilter(inputs=inputs) field_interpolation_filter.UpdateAttributes() document.AllAlgorithms.Add(field_interpolation_filter)Then you will have the same grid and array ordering, and can overlay the T1w image and the field (or display them side by side with identical aspect ratio).
Note if your simulation grid is not uniform, matplotlib will depict the anatomy in a distorted way.
Btw, you could also interpolate in the other direction ... -
You don't have the same resolution (grid spacing) in the T1 image and the efield...!
To visualize the two "fields" on the same grid, drag the T1w image to the analysis explorer, and interpolate it on the e-field grid.The interpolator has two inputs:
- T1w image (input field)
- E-field (target grid)
Here is a snippet to show this could be done in Python
import s4l_v1.analysis as analysis import s4l_v1.document as document import s4l_v1.model as model # Add a new ModelToGridFilter inputs = [] model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs) model_to_grid_filter.Name = "Image" model_to_grid_filter.Entity = model.AllEntities()["T1w-image"] # The model entity corresponding to the T1w-image model_to_grid_filter.UpdateAttributes() # Add a new SimulationExtractor simulation = document.AllSimulations["EM"] simulation_extractor = simulation.Results() # Add a new EmSensorExtractor em_sensor_extractor = simulation_extractor["Overall Field"] em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All" # Add a new FieldInterpolationFilter inputs = [model_to_grid_filter.Outputs[""], em_sensor_extractor.Outputs["EM E(x,y,z,f0)"]] field_interpolation_filter = analysis.core.FieldInterpolationFilter(inputs=inputs) field_interpolation_filter.UpdateAttributes() document.AllAlgorithms.Add(field_interpolation_filter)Then you will have the same grid and array ordering, and can overlay the T1w image and the field (or display them side by side with identical aspect ratio).
Note if your simulation grid is not uniform, matplotlib will depict the anatomy in a distorted way.
Btw, you could also interpolate in the other direction ... -
what do you mean by "The T1 and the electric field seem to be in opposite directions along the x-axis"?
- the interpolated T1 image and the E-field?
- is the interpolated T1 image NOT shown in the same position/orientation as the original T1 image from the modeler?
- is the original T1 image (which you imported in the modeler) aligned/shown in the same position/orientation as the model?
It would help if you would add screenshots to your question, with explanations of what each screenshot is trying to demonstrate.
-
what do you mean by "The T1 and the electric field seem to be in opposite directions along the x-axis"?
- the interpolated T1 image and the E-field?
- is the interpolated T1 image NOT shown in the same position/orientation as the original T1 image from the modeler?
- is the original T1 image (which you imported in the modeler) aligned/shown in the same position/orientation as the model?
It would help if you would add screenshots to your question, with explanations of what each screenshot is trying to demonstrate.
@bryn What I mean is that the slice with index 0 in the T1 image and the slice with index 1 in the electric field are not on the same side of the head model along the x-axis. Additionally, I would like to know if it is possible to export the interpolated T1 image. I need the T1 structural image that corresponds to the simulated electric field to create a dataset.
-
but are you using a slice field viewer to show the interpolated MRI slice? If not you basically are looking at the original MRI, with different grid and data ordering.
-
I sent you code (it can also be done in GUI) above to interpolate the T1w image onto the E-field/simulation grid. If you do that you can export it in the same format as you export the E-field (not as a nii.gz file, since the simulation grid is rectilinear which is not supported by nifti files).
-
@bryn Is the following code correct? I'm not sure if the input is correct.
inputs1 = [model_to_grid_filter.Outputs[""]] field_data_text_exporter1 = analysis.exporters.FieldDataTextExporter(inputs=inputs1) field_data_text_exporter1.FileName = output_file1 field_data_text_exporter1.UpdateAttributes() document.AllAlgorithms.Add(field_data_text_exporter1) field_data_text_exporter1.Update(overwrite=True) -
this does not include the interpolation...
to access an output of an analysis algorithm you can use named outputs or indexed outputs, in your snippet above it may be easier to use the index '0'.
inputs1 = [model_to_grid_filter.Outputs[0]] -
what does your code look like? I create the pipeline in the GUI and it seems to run fine for me (v9.0.0). The "To Python" function yields
# This script was auto-generated by Sim4Life version 9.0.1 import numpy import s4l_v1.analysis as analysis import s4l_v1.document as document import s4l_v1.model as model import s4l_v1.units as units from s4l_v1 import ReleaseVersion from s4l_v1 import Unit try: # Define the version to use for default values ReleaseVersion.set_active(ReleaseVersion.version9_0) # Creating the analysis pipeline # Adding a new ModelToGridFilter inputs = [] model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs) model_to_grid_filter.Name = "IXI025-Guys-0852-T1" model_to_grid_filter.Entity = model.AllEntities()["IXI025-Guys-0852-T1"] model_to_grid_filter.UpdateAttributes() document.AllAlgorithms.Add(model_to_grid_filter) # Adding a new SimulationExtractor simulation = document.AllSimulations["EM"] simulation_extractor = simulation.Results() # Adding a new EmSensorExtractor em_sensor_extractor = simulation_extractor["Overall Field"] em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All" document.AllAlgorithms.Add(em_sensor_extractor) # Adding a new FieldInterpolationFilter inputs = [model_to_grid_filter.Outputs[0], em_sensor_extractor.Outputs["EM E(x,y,z,f0)"]] field_interpolation_filter = analysis.core.FieldInterpolationFilter(inputs=inputs) field_interpolation_filter.UpdateAttributes() document.AllAlgorithms.Add(field_interpolation_filter) # Adding a new FieldDataTextExporter inputs = [field_interpolation_filter.Outputs[""]] field_data_text_exporter = analysis.exporters.FieldDataTextExporter(inputs=inputs) field_data_text_exporter.FileName = u"D:\\temp\\ExportedFieldData-T1w.txt" field_data_text_exporter.UpdateAttributes() document.AllAlgorithms.Add(field_data_text_exporter) except Exception as exc: import traceback traceback.print_exc() # Reset active version to default ReleaseVersion.reset() raise(exc) -
@bryn This is my code, I use v8.2.2
def run_simulation(value1, value2): print(f"Running simulation for values: {value1}_{value2}") # Define the version to use for default values # s2=time.time() #simulation # Creating the simulation simulation1 = emlf.ElectroQsOhmicSimulation() simulation1.Name =f"{value1}_{value2}" # Mapping the components and entities component__plane_x = simulation1.AllComponents["Plane X+"] component__plane_x = simulation1.AllComponents["Plane X-"] component__background = simulation1.AllComponents["Background"] component__plane_y = simulation1.AllComponents["Plane Y+"] component__plane_y = simulation1.AllComponents["Plane Y-"] component__plane_z = simulation1.AllComponents["Plane Z+"] component__plane_z = simulation1.AllComponents["Plane Z-"] component__overall_field = simulation1.AllComponents["Overall Field"] entity_0 = model.AllEntities()["Cz_Cylinder 1"] entity_1 = model.AllEntities()["Fpz_Cylinder 1"] entity_2 = model.AllEntities()["Fp1_Cylinder 1"] entity_3 = model.AllEntities()["AF7_Cylinder 1"] entity_4 = model.AllEntities()["F7_Cylinder 1"] entity_5 = model.AllEntities()["FT7_Cylinder 1"] entity_6 = model.AllEntities()["T7_Cylinder 1"] entity_7 = model.AllEntities()["TP7_Cylinder 1"] entity_8 = model.AllEntities()["P7_Cylinder 1"] entity_9 = model.AllEntities()["PO7_Cylinder 1"] entity_10 = model.AllEntities()["O1_Cylinder 1"] entity_11 = model.AllEntities()["Oz_Cylinder 1"] entity_12 = model.AllEntities()["O2_Cylinder 1"] entity_13 = model.AllEntities()["PO8_Cylinder 1"] entity_14 = model.AllEntities()["P8_Cylinder 1"] entity_15 = model.AllEntities()["TP8_Cylinder 1"] entity_16 = model.AllEntities()["T8_Cylinder 1"] entity_17 = model.AllEntities()["FT8_Cylinder 1"] entity_18 = model.AllEntities()["F8_Cylinder 1"] entity_19 = model.AllEntities()["AF8_Cylinder 1"] entity_20 = model.AllEntities()["Fp2_Cylinder 1"] entity_21 = model.AllEntities()["AFz_Cylinder 1"] entity_22 = model.AllEntities()["AF3_Cylinder 1"] entity_23 = model.AllEntities()["F5_Cylinder 1"] entity_24 = model.AllEntities()["FC5_Cylinder 1"] entity_25 = model.AllEntities()["C5_Cylinder 1"] entity_26 = model.AllEntities()["CP5_Cylinder 1"] entity_27 = model.AllEntities()["P5_Cylinder 1"] entity_28 = model.AllEntities()["PO3_Cylinder 1"] entity_29 = model.AllEntities()["POz_Cylinder 1"] entity_30 = model.AllEntities()["PO4_Cylinder 1"] entity_31 = model.AllEntities()["P6_Cylinder 1"] entity_32 = model.AllEntities()["CP6_Cylinder 1"] entity_33= model.AllEntities()["C6_Cylinder 1"] entity_34 = model.AllEntities()["FC6_Cylinder 1"] entity_35 = model.AllEntities()["F6_Cylinder 1"] entity_36 = model.AllEntities()["AF4_Cylinder 1"] entity_37 = model.AllEntities()["Fz_Cylinder 1"] entity_38 = model.AllEntities()["F1_Cylinder 1"] entity_39 = model.AllEntities()["F3_Cylinder 1"] entity_40 = model.AllEntities()["FC3_Cylinder 1"] entity_41 = model.AllEntities()["C3_Cylinder 1"] entity_42 = model.AllEntities()["CP3_Cylinder 1"] entity_43 = model.AllEntities()["P3_Cylinder 1"] entity_44 = model.AllEntities()["P1_Cylinder 1"] entity_45 = model.AllEntities()["Pz_Cylinder 1"] entity_46 = model.AllEntities()["P2_Cylinder 1"] entity_47 = model.AllEntities()["P4_Cylinder 1"] entity_48 = model.AllEntities()["CP4_Cylinder 1"] entity_49 = model.AllEntities()["C4_Cylinder 1"] entity_50 = model.AllEntities()["FC4_Cylinder 1"] entity_51 = model.AllEntities()["F4_Cylinder 1"] entity_52 = model.AllEntities()["F2_Cylinder 1"] entity_53 = model.AllEntities()["FCz_Cylinder 1"] entity_54 = model.AllEntities()["FC1_Cylinder 1"] entity_55 = model.AllEntities()["C1_Cylinder 1"] entity_56 = model.AllEntities()["CP1_Cylinder 1"] entity_57 = model.AllEntities()["CPz_Cylinder 1"] entity_58 = model.AllEntities()["CP2_Cylinder 1"] entity_59 = model.AllEntities()["C2_Cylinder 1"] entity_60 = model.AllEntities()["FC2_Cylinder 1"] entity__bone_cortical = model.AllEntities()["Bone_cortical"] entity__vein = model.AllEntities()["Vein"] entity__other_tissues = model.AllEntities()["Other_tissues"] entity__nerve_cranial_ii_optic = model.AllEntities()["Nerve_cranial_II_optic"] entity__muscle_ocular = model.AllEntities()["Muscle_ocular"] entity__bone_cancellous = model.AllEntities()["Bone_cancellous"] entity__spinal_cord = model.AllEntities()["Spinal_cord"] entity__skin = model.AllEntities()["Skin"] entity__cerebrum_white_matter = model.AllEntities()["Cerebrum_white_matter"] entity__eyes = model.AllEntities()["Eyes"] entity__air_internal = model.AllEntities()["Air_internal"] entity__mucosa = model.AllEntities()["Mucosa"] entity__cerebrospinal_fluid = model.AllEntities()["Cerebrospinal_fluid"] entity__dura = model.AllEntities()["Dura"] entity__artery = model.AllEntities()["Artery"] entity__cerebrum_grey_matter = model.AllEntities()["Cerebrum_grey_matter"] # 创建一个映射表，将数字索引与实体名称关联起来 entity_mapping = { 0: "Cz_Cylinder 1", 1: "Fpz_Cylinder 1", 2: "Fp1_Cylinder 1", 3: "AF7_Cylinder 1", 4: "F7_Cylinder 1", 5: "FT7_Cylinder 1", 6: "T7_Cylinder 1", 7: "TP7_Cylinder 1", 8: "P7_Cylinder 1", 9: "PO7_Cylinder 1", 10: "O1_Cylinder 1", 11: "Oz_Cylinder 1", 12: "O2_Cylinder 1", 13: "PO8_Cylinder 1", 14: "P8_Cylinder 1", 15: "TP8_Cylinder 1", 16: "T8_Cylinder 1", 17: "FT8_Cylinder 1", 18: "F8_Cylinder 1", 19: "AF8_Cylinder 1", 20: "Fp2_Cylinder 1", 21: "AFz_Cylinder 1", 22: "AF3_Cylinder 1", 23: "F5_Cylinder 1", 24: "FC5_Cylinder 1", 25: "C5_Cylinder 1", 26: "CP5_Cylinder 1", 27: "P5_Cylinder 1", 28: "PO3_Cylinder 1", 29: "POz_Cylinder 1", 30: "PO4_Cylinder 1", 31: "P6_Cylinder 1", 32: "CP6_Cylinder 1", 33: "C6_Cylinder 1", 34: "FC6_Cylinder 1", 35: "F6_Cylinder 1", 36: "AF4_Cylinder 1", 37: "Fz_Cylinder 1", 38: "F1_Cylinder 1", 39: "F3_Cylinder 1", 40: "FC3_Cylinder 1", 41: "C3_Cylinder 1", 42: "CP3_Cylinder 1", 43: "P3_Cylinder 1", 44: "P1_Cylinder 1", 45: "Pz_Cylinder 1", 46: "P2_Cylinder 1", 47: "P4_Cylinder 1", 48: "CP4_Cylinder 1", 49: "C4_Cylinder 1", 50: "FC4_Cylinder 1", 51: "F4_Cylinder 1", 52: "F2_Cylinder 1", 53: "FCz_Cylinder 1", 54: "FC1_Cylinder 1", 55: "C1_Cylinder 1", 56: "CP1_Cylinder 1", 57: "CPz_Cylinder 1", 58: "CP2_Cylinder 1", 59: "C2_Cylinder 1", 60: "FC2_Cylinder 1" } # 确保输入的值是整数类型 value1 = int(value1) value2 = int(value2) # 根据输入的数字索引获取对应的实体名称 entity_name1 = entity_mapping.get(value1) entity_name2 = entity_mapping.get(value2) # 检查实体名称是否存在 if entity_name1 is None: raise KeyError(f"Entity for value {value1} does not exist") if entity_name2 is None: raise KeyError(f"Entity for value {value2} does not exist") # 获取对应的实体对象 entity_value1 = model.AllEntities()[entity_name1] entity_value2 = model.AllEntities()[entity_name2] # 打印结果，确认是否正确 print(f"Entity for value1 ({value1}): {entity_value1}") print(f"Entity for value2 ({value2}): {entity_value2}") # Editing QuasiStaticSetupSettings "Setup quasi_static_setup_settings = [x for x in simulation1.AllSettings if isinstance(x, emlf.QuasiStaticSetupSettings) and x.Name == "Setup"][0] quasi_static_setup_settings.Frequency = 2000.0, units.Hz # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__skin] mat = database["IT'IS LF 4.2"]["Skin"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Skin" material_settings.MassDensity = 1109.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.1482971014492752, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 1135.619382618975 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__air_internal] mat = database["IT'IS LF 4.2"]["Air"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Air 1" material_settings.MassDensity = 1.164091552938177, Unit("kg/m^3") material_settings.Name = "Air 1" simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__cerebrospinal_fluid] mat = database["IT'IS LF 4.2"]["Cerebrospinal Fluid"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Cerebrospinal Fluid" material_settings.MassDensity = 1007.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 1.878999709695023, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 108.9999972649334 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__cerebrum_white_matter] mat = database["IT'IS LF 4.2"]["Brain (White Matter)"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Brain (White Matter)" material_settings.MassDensity = 1041.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.3479543931346832, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 69810.68883114036 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__mucosa] mat = database["IT'IS LF 4.2"]["Mucous Membrane"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Mucous Membrane" material_settings.MassDensity = 1102.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.4610075264456888, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 434932.19242741907 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__dura] mat = database["IT'IS LF 4.2"]["Dura"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Dura" material_settings.MassDensity = 1174.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.06, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 5343.990898234167 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__eyes] mat = database["IT'IS LF 4.2"]["Eye (Aqueous Humor)"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Eye (Aqueous Humor)" material_settings.MassDensity = 993.7770167788401, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 1.878999709695023, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 108.9999972649334 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__artery, entity__vein] mat = database["IT'IS LF 4.2"]["Blood"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Blood" material_settings.MassDensity = 1049.75, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.6624597361833767, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 5258.608390020375 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__spinal_cord] mat = database["IT'IS LF 4.2"]["Spinal Cord"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Spinal Cord" material_settings.MassDensity = 1075.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.6109538492063492, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 69911.4914652573 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__nerve_cranial_ii_optic] mat = database["IT'IS LF 4.2"]["Nerve"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Nerve" material_settings.MassDensity = 1075.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.3479543931346832, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 69911.4914652573 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__muscle_ocular] mat = database["IT'IS LF 4.2"]["Muscle"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Muscle" material_settings.MassDensity = 1090.4, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.4610075264456888, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 434932.19242741907 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__bone_cortical] mat = database["IT'IS LF 4.2"]["Bone (Cortical)"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Bone (Cortical)" material_settings.MassDensity = 1908.0, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.00630199709513435, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 2702.3711256306647 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__cerebrum_grey_matter] mat = database["IT'IS LF 4.2"]["Brain (Grey Matter)"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Brain (Grey Matter)" material_settings.MassDensity = 1044.5, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.4190548817650446, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 164062.99316639948 simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__other_tissues] material_settings.ElectricProps.Conductivity = 0.087, Unit("S/m") simulation1.Add(material_settings, components) # Adding a new MaterialSettings material_settings = emlf.MaterialSettings() components = [entity__bone_cancellous] mat = database["IT'IS LF 4.2"]["Bone (Cancellous)"] if mat is not None: simulation1.LinkMaterialWithDatabase(material_settings, mat) else: # Fallback if material is not found material_settings.Name = "Bone (Cancellous)" material_settings.MassDensity = 1178.333333333333, Unit("kg/m^3") material_settings.ElectricProps.Conductivity = 0.08045952772338552, Unit("S/m") material_settings.ElectricProps.RelativePermittivity = 12320.035797440474 simulation1.Add(material_settings, components) # Editing BoundarySettings "Boundary Settings boundary_settings = [x for x in simulation1.AllSettings if isinstance(x, emlf.BoundarySettings) and x.Name == "Boundary Settings"][0] components = [component__plane_x, component__plane_x, component__plane_y, component__plane_y, component__plane_z, component__plane_z] simulation1.Add(boundary_settings, components) boundary_settings.BoundaryType = boundary_settings.BoundaryType.enum.Flux # Adding a new BoundarySettings boundary_settings = emlf.BoundarySettings() components = [entity_value1] boundary_settings.Name = "Boundary Settings 1" boundary_settings.DirichletValue = 1.0, units.Volts simulation1.Add(boundary_settings, components) # Adding a new BoundarySettings boundary_settings = emlf.BoundarySettings() components = [entity_value2] boundary_settings.Name = "Boundary Settings 2" boundary_settings.DirichletValue = -1.0, units.Volts simulation1.Add(boundary_settings, components) # Editing GlobalGridSettings "Grid (Empty) global_grid_settings = simulation1.GlobalGridSettings global_grid_settings.PaddingMode = global_grid_settings.PaddingMode.enum.Manual global_grid_settings.BottomPadding = numpy.array([0.0, 0.0, 0.0]), units.MilliMeters global_grid_settings.TopPadding = numpy.array([0.0, 0.0, 0.0]), units.MilliMeters # Adding a new ManualGridSettings manual_grid_settings = simulation1.AddManualGridSettings([entity_0, entity_1, entity_2,entity_3,entity_4,entity_5,entity_6,entity_7,entity_8,entity_9,entity_10, entity_11,entity_12,entity_13,entity_14,entity_15,entity_16,entity_17,entity_18,entity_19,entity_20, entity_21,entity_22,entity_23,entity_24,entity_25,entity_26,entity_27,entity_28,entity_29,entity_30, entity_31,entity_32,entity_33,entity_34,entity_35,entity_36,entity_37,entity_38,entity_39,entity_40, entity_41,entity_42,entity_43,entity_44,entity_45,entity_46,entity_47,entity_48,entity_49,entity_50, entity_51,entity_52,entity_53,entity_54,entity_55,entity_56,entity_57,entity_58,entity_59,entity_60]) manual_grid_settings.MaxStep = numpy.array([1.0, 1.0, 1.0]), units.MilliMeters manual_grid_settings.Resolution = numpy.array([0.3, 0.3, 0.3]), units.MilliMeters manual_grid_settings.Priority = 0.0 # Adding a new ManualGridSettings manual_grid_settings = simulation1.AddManualGridSettings([entity__air_internal, entity__artery, entity__bone_cancellous, entity__bone_cortical, entity__cerebrospinal_fluid, entity__cerebrum_grey_matter, entity__cerebrum_white_matter, entity__dura, entity__eyes, entity__mucosa, entity__muscle_ocular, entity__nerve_cranial_ii_optic, entity__other_tissues, entity__skin, entity__spinal_cord, entity__vein]) manual_grid_settings.MaxStep = numpy.array([1.0, 1.0, 1.0]), units.MilliMeters manual_grid_settings.Resolution = numpy.array([0.3, 0.3, 0.3]), units.MilliMeters manual_grid_settings.Priority = 0.0 # Editing AutomaticVoxelerSettings "Automatic Voxeler Settings automatic_voxeler_settings = [x for x in simulation1.AllSettings if isinstance(x, emlf.AutomaticVoxelerSettings) and x.Name == "Automatic Voxeler Settings"][0] components = [entity__air_internal, entity__artery, entity__bone_cancellous, entity__bone_cortical, entity__cerebrospinal_fluid, entity__cerebrum_grey_matter, entity__cerebrum_white_matter, entity__dura, entity__eyes, entity__mucosa, entity__muscle_ocular, entity__nerve_cranial_ii_optic, entity__other_tissues, entity__skin, entity__spinal_cord, entity__vein] simulation1.Add(automatic_voxeler_settings, components) automatic_voxeler_settings.Priority = 1 # Adding a new AutomaticVoxelerSettings automatic_voxeler_settings = emlf.AutomaticVoxelerSettings() #components = [entity_f5__cylinder1, entity_p5__cylinder1] components = [entity_value1, entity_value2] automatic_voxeler_settings.Name = "electrodes" simulation1.Add(automatic_voxeler_settings, components) automatic_voxeler_settings.Priority = 0 # Editing SolverSettings "Solver solver_settings = simulation1.SolverSettings solver_settings.NumberOfProcesses = 2 # Update the materials with the new frequency parameters simulation1.UpdateAllMaterials() # Update the grid with the new parameters simulation1.UpdateGrid() # simulation1.CreateVoxels(r"F:\Huilin\sim4life_pro\script\head_TI\head_TI-4.smash") # Add the simulation to the UI document.AllSimulations.Add( simulation1 ) simulation1.RunSimulation(wait=True) # Output directory output_dir = r"G:\MHL\TI_head\IXI025\interprate\tacs_2000Hz" # Output file output_file1 = os.path.join(output_dir, f"{value1}_{value2}.txt") output_file2 = os.path.join(output_dir, f"{value1}_{value2}_T1_interpolation.txt") # Add a new ModelToGridFilter inputs = [] model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs) model_to_grid_filter.Name = "Image" model_to_grid_filter.Entity = model.AllEntities()["IXI025-Guys-0852-T1"] # The model entity corresponding to the T1w-image model_to_grid_filter.UpdateAttributes() simulation1 = document.AllSimulations[f"{value1}_{value2}"] simulation_extractor1 = simulation1.Results() # Adding a new EmSensorExtractor em_sensor_extractor1 = simulation_extractor1["Overall Field"] em_sensor_extractor1.FrequencySettings.ExtractedFrequency = u"All" document.AllAlgorithms.Add(em_sensor_extractor1) # Adding a new CurrentExtractor inputs1 = [em_sensor_extractor1.Outputs["EM Potential(x,y,z,f0)"], em_sensor_extractor1.Outputs["J(x,y,z,f0)"]] current_extractor1 = analysis.extractors.CurrentExtractor(inputs=inputs1) current_extractor1.IsoSurfacePerCentThreshold = 30 current_extractor1.UpdateAttributes() document.AllAlgorithms.Add(current_extractor1) # Adding a new DataTableHTMLViewer output1 = current_extractor1.Outputs["Total Flux(J(x,y,z,f0))"] output1.Update() comp1 = output1.Data.GetComponent(0) flux1 = np.real(comp1)[0] inputs1 = [em_sensor_extractor1.Outputs["EM E(x,y,z,f0)"]] user_defined_field_normalizer1 = analysis.field.UserDefinedFieldNormalizer(inputs=inputs1) user_defined_field_normalizer1.Target.Value = 1.0 / (flux1 * 1000) user_defined_field_normalizer1.Name = f"{value1}_{value2}" user_defined_field_normalizer1.UpdateAttributes() document.AllAlgorithms.Add(user_defined_field_normalizer1) # Adding a new FieldDataTextExporter inputs1 = [user_defined_field_normalizer1.Outputs["EM E(x,y,z,f0)"]] field_data_text_exporter1 = analysis.exporters.FieldDataTextExporter(inputs=inputs1) field_data_text_exporter1.FileName = output_file1 field_data_text_exporter1.UpdateAttributes() document.AllAlgorithms.Add(field_data_text_exporter1) field_data_text_exporter1.Update(overwrite=True) inputs2 = [model_to_grid_filter.Outputs[0]] field_data_text_exporter2 = analysis.exporters.FieldDataTextExporter(inputs=inputs2) field_data_text_exporter2.FileName = output_file2 field_data_text_exporter2.UpdateAttributes() document.AllAlgorithms.Add(field_data_text_exporter2) field_data_text_exporter2.Update(overwrite=True)
