Merge pull request #49 from SFB-ELAINE/fix_edge_refinement

Fix edge refinement

.github/ISSUE_TEMPLATE

@@ -0,0 +1,16 @@
+* OSS-DBS version (either Lead-DBS or development, then please add which branch / commit):
+* If installed through Lead-DBS, please state the version:
+* Python version:
+* Operating System:
+
+### Description
+
+Describe what you were trying to get done.
+Tell us what happened, what went wrong, and what you expected to happen.
+
+### What I Did
+
+```
+Paste the command(s) you ran and the output.
+If there was a crash, please include the traceback here.
+```

.github/workflows/ci.yml

@@ -27,12 +27,8 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
         platform: [ubuntu-latest, macos-latest]
-        exclude:
-          - platform: macos-latest
-            python-version: 3.8
-
     steps:
       - uses: actions/checkout@v4
 

.gitignore

@@ -28,4 +28,11 @@ impedance.csv
 *.nii.gz
 *.nii
 success_*.txt
-fail_*.txt
+fail_*.txt
+test_input_cases_results.csv
+*.stp
+ngmesh.ini
+examples/StimulationSignalAPI/*.pdf
+x86_64/
+
+meshsizes.txt

.pre-commit-config.yaml

@@ -13,13 +13,9 @@ repos:
     hooks:
       - id: validate-pyproject
 
-  - repo: https://github.com/charliermarsh/ruff-pre-commit
-    rev: v0.0.282
+  - repo: https://github.com/astral-sh/ruff-pre-commit
+    rev: v0.6.2
     hooks:
       - id: ruff
         args: [--fix]
-
-  - repo: https://github.com/psf/black
-    rev: 23.7.0
-    hooks:
-      - id: black
+      - id: ruff-format

examples/Collapse_VTA/inputCollapse.json

@@ -136,6 +136,6 @@
   },
   "OutputPath": "Results",
   "SaveImpedance": true,
-  "ActivationThresholdVTA": 0.2,
+  "ActivationThresholdVTA[V-per-m]": 0.2,
   "ExportVTK": true
 }

examples/ConvergenceStudy/FiguresReproducibility/MeshRefinement/generate_meshes.py

@@ -0,0 +1,151 @@
+# gold standard:
+# impedance error 0.1% and maximum element size equal
+# to mesh size
+
+import json
+import logging
+import os
+from copy import deepcopy
+
+import ngsolve
+import numpy as np
+
+import ossdbs
+
+
+def remove_file_handler(logger):
+    """Remove file handler from logger instance."""
+    for h in logger.handlers:
+        if isinstance(h, logging.FileHandler):
+            logger.removeHandler(h)
+
+
+def save_input_dict(base_input_dict):
+    """Save input dictionary for convergence run."""
+    input_dict = deepcopy(base_input_dict)
+    input_dict["OutputPath"] = "./"
+
+    with open(
+        os.path.join(base_input_dict["OutputPath"], "input_dict.json"), "w"
+    ) as fp:
+        json.dump(input_dict, fp, indent=2)
+
+
+BASE_CONTACT = {
+    "Contact_ID": 1,
+    "Active": False,
+    "Current[A]": 0.0,
+    "Voltage[V]": 0.0,
+    "Floating": False,
+}
+
+ossdbs.set_logger(logging.INFO)
+_logger = logging.getLogger("ossdbs")
+
+electrode_name = "Medtronic3389"
+
+with open("oss_dbs_parameters.json") as fp:
+    base_input_dict = json.load(fp)
+
+
+# add electrode
+base_input_dict["Electrodes"][0]["Name"] = electrode_name
+
+# initially no adaptive refinement
+base_input_dict["Mesh"]["AdaptiveMeshRefinement"] = {}
+base_input_dict["Mesh"]["AdaptiveMeshRefinement"]["Active"] = False
+
+# mri image
+mri_image, _ = ossdbs.load_images(base_input_dict)
+# pathway info
+pw = ossdbs.point_analysis.Pathway(base_input_dict["PointModel"]["Pathway"]["FileName"])
+pw.write_netgen_meshsize_file(
+    meshsize=min(mri_image.voxel_sizes), filename="meshsizes.txt"
+)
+
+# first refinement level: Default
+base_input_dict["Mesh"]["MeshingHypothesis"]["Type"] = "Default"
+base_input_dict["Mesh"]["MeshingHypothesis"]["MaxMeshSize"] = 1e6
+base_input_dict["OutputPath"] = "Results_PAM_default"
+electrodes = ossdbs.api.generate_electrodes(base_input_dict)
+brain_model = ossdbs.api.build_brain_model(base_input_dict, mri_image)
+try:
+    geometry = ossdbs.ModelGeometry(brain_model, electrodes)
+except RuntimeError:
+    brain_model = ossdbs.api.build_brain_model(
+        base_input_dict, mri_image, rotate_initial_geo=True
+    )
+    geometry = ossdbs.ModelGeometry(brain_model, electrodes)
+
+
+solver = ossdbs.api.prepare_solver(base_input_dict)
+dielectric_properties = ossdbs.api.prepare_dielectric_properties(base_input_dict)
+materials = base_input_dict["MaterialDistribution"]["MRIMapping"]
+conductivity = ossdbs.ConductivityCF(
+    mri_image,
+    brain_model.brain_region,
+    dielectric_properties,
+    materials,
+    geometry.encapsulation_layers,
+    complex_data=base_input_dict["EQSMode"],
+)
+
+
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+
+# refinement level: default and refinement around pathway
+base_input_dict["Mesh"]["MeshingHypothesis"]["MeshSizeFilename"] = "meshsizes.txt"
+base_input_dict["OutputPath"] = "Results_PAM_default_meshsize"
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+# reset mesh size file name
+base_input_dict["Mesh"]["MeshingHypothesis"]["MeshSizeFilename"] = ""
+
+
+# fourth refinement: material refinement
+base_input_dict["OutputPath"] = "Results_PAM_material_refinement"
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    # two refinement steps
+    volume_conductor.refine_mesh_by_material(2)
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+
+# sixth refinement: more edge refinement
+base_input_dict["Mesh"]["MeshingHypothesis"]["Type"] = "Default"
+lead_diameter = ossdbs.electrodes.default_electrode_parameters[
+    electrode_name
+].lead_diameter
+perimeter = np.pi * lead_diameter
+edge_size = perimeter / 50.0
+print(edge_size)
+base_input_dict["Electrodes"][0]["Contacts"][0]["MaxMeshSizeEdge"] = edge_size
+base_input_dict["Electrodes"][0]["Contacts"][1]["MaxMeshSizeEdge"] = edge_size
+base_input_dict["OutputPath"] = "Results_PAM_fine_edge_refinement"
+ossdbs.api.set_contact_and_encapsulation_layer_properties(base_input_dict, geometry)
+
+
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    geo = volume_conductor.mesh.ngsolvemesh.ngmesh.GetGeometry()
+    for edge in geo.shape.edges:
+        print(edge.name, edge.maxh)
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()

examples/ConvergenceStudy/FiguresReproducibility/MeshRefinement/oss_dbs_parameters.json

@@ -0,0 +1,150 @@
+{
+    "BrainRegion": {
+        "Center": {"x[mm]": 13.5035, "y[mm]": -9.49908, "z[mm]": -9.87216},
+        "Dimension": {"x[mm]": 60, "y[mm]": 60, "z[mm]": 80},
+        "Shape": "Ellipsoid"
+    },
+    "Electrodes": [
+        {
+            "Name": "",
+            "Rotation[Degrees]": 0.15176099498051523,
+            "Direction": {
+                "x[mm]": 0.0,
+                "y[mm]": 0.0,
+                "z[mm]": 1.0
+            },
+            "TipPosition": {
+                "x[mm]": 12.956271802141353,
+                "y[mm]": -9.901870551007098,
+                "z[mm]": -12.01710780157648
+            },
+            "EncapsulationLayer": {
+                "Thickness[mm]": 0.0,
+                "Material": "Gray matter",
+                "DielectricModel": "ColeCole4",
+                "DielectricParameters": null,
+                "MaxMeshSize": 0.1
+            },
+            "Contacts": [
+                {
+                    "Contact_ID": 1,
+                    "Active": false,
+                    "Current[A]": 0.0,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 2,
+                    "Active": false,
+                    "Current[A]": 0.012,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 3,
+                    "Active": false,
+                    "Current[A]": -0.004,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 4,
+                    "Active": false,
+                    "Current[A]": -0.008,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                }
+            ]
+        }
+    ],
+    "Surfaces": [
+        {"Name": "BrainSurface", "Active": true, "Current[A]": 0.0, "Voltage[V]": 0.0}
+    ],
+    "MaterialDistribution": {
+        "MRIPath": "../../PAM/segmask.nii.gz",
+        "MRIMapping": {
+            "Blood": 4,
+            "Gray matter": 1,
+            "White matter": 2,
+            "CSF": 3,
+            "Unknown": 0
+        },
+        "DiffusionTensorActive": false,
+        "DTIPath": ""
+    },
+    "DielectricModel": {"Type": "ColeCole4", "CustomParameters": null},
+    "Mesh": {
+        "LoadMesh": false,
+        "LoadPath": "",
+        "MeshingHypothesis": {"Type": "Default"},
+        "HPRefinement": {"Active": false, "Levels": 2, "Factor": 0.125},
+        "AdaptiveMeshRefinement": {
+            "Active": false,
+            "MaxIterations": 10,
+            "ErrorTolerance": 0.1
+        },
+        "MeshSize": {"Edges": {}, "Faces": {}, "Volumes": {}},
+        "SaveMesh": false,
+        "SavePath": "mesh"
+    },
+    "EQSMode": false,
+    "FEMOrder": 2,
+    "StimulationSignal": {
+        "Type": "Rectangle",
+        "ListOfFrequencies": [10000],
+        "Frequency[Hz]": 130.0,
+        "PulseWidth[us]": 30.0,
+        "PulseTopWidth[us]": 0.0,
+        "CounterPulseWidth[us]": 0.0,
+        "InterPulseWidth[us]": 0.0,
+        "SpectrumMode": "OctaveBand",
+        "CounterAmplitude": 1.0,
+        "CutoffFrequency": 250000.0,
+        "CurrentControlled": true
+    },
+    "Solver": {
+        "Type": "CG",
+        "Preconditioner": "local",
+        "PreconditionerKwargs": {},
+        "MaximumSteps": 5000,
+        "Precision": 1e-10
+    },
+    "PointModel": {
+        "Pathway": {
+            "Active": true,
+            "FileName": "../../PAM/Allocated_axons.h5",
+            "ExportField": false
+        },
+        "Lattice": {
+            "Active": false,
+            "Center": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
+            "Shape": {"x": 10, "y": 10, "z": 10},
+            "Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
+            "PointDistance[mm]": 0.1,
+            "CollapseVTA": false,
+            "ExportField": false
+        },
+        "VoxelLattice": {
+            "Active": false,
+            "Shape": {"x": 10, "y": 10, "z": 10},
+            "TimeDomain": false,
+            "ExportField": false
+        }
+    },
+    "OutputPath": "Results_rh",
+    "ComputeImpedance": false,
+    "ExportVTK": true,
+    "ExportFrequency": null,
+    "ExportElectrode": false,
+    "ModelSide": 0,
+    "CalcAxonActivation": true,
+    "ActivationThresholdVTA[V-per-m]": 200.0,
+    "FailFlag": "rh",
+    "OutOfCore": false,
+    "PathwayFile": "../../PAM/Allocated_axons_parameters.json",
+    "StimSets": {"Active": false, "StimSetsFile": null}
+}