diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 7a73a41b..00000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,2 +0,0 @@ -{ -} \ No newline at end of file diff --git a/Makefile b/Makefile index 41778a41..5910cdca 100644 --- a/Makefile +++ b/Makefile @@ -4,11 +4,11 @@ help: @echo "Available commands are: \n*lint, lint-check, spellcheck, test, test-unit, test-integration docs" lint: - poetry run isort src tests - poetry run black src tests - poetry run flake8 src tests - poetry run mypy src - poetry run codespell src + poetry run isort src tests docs/examples + poetry run black src tests docs/examples + poetry run flake8 src tests docs/examples + poetry run mypy src docs/examples + poetry run codespell src docs/examples poetry run pydocstyle src poetry run pyright @@ -16,11 +16,11 @@ spellcheck: poetry run pylint --disable all --enable spelling --spelling-dict en_US --spelling-private-dict-file=whitelist.txt src lint-check: - poetry run isort --check src tests - poetry run black --check src tests - poetry run flake8 src tests - poetry run mypy src - poetry run codespell src + poetry run isort --check src tests docs/examples + poetry run black --check src tests docs/examples + poetry run flake8 src tests docs/examples + poetry run mypy src docs/examples + poetry run codespell src docs/examples poetry run pydocstyle src poetry run pyright --warnings diff --git a/docs/examples/plot_3d.py b/docs/examples/plot_3d.py index 5289d2b1..a90bba62 100644 --- a/docs/examples/plot_3d.py +++ b/docs/examples/plot_3d.py @@ -3,7 +3,8 @@ Visualizing 3D results with Napari ==================================== -This example demonstrates how to render a steady state result for a 3D scenario using [napari](https://napari.org/stable/). +This example demonstrates how to render a steady state result for +a 3D scenario using [napari](https://napari.org/stable/). Running such simulations is computationally expensive and can take a long time to complete. For this reason, we recommend running this simulation on an external machine, @@ -15,13 +16,12 @@ # %% # The following step downloads and loads a simulation executed on an external machine. import pooch + import neurotechdevkit as ndk -URL = 'https://neurotechdevkit.s3.us-west-2.amazonaws.com/result-scenario-2-3d.tz' +URL = "https://neurotechdevkit.s3.us-west-2.amazonaws.com/result-scenario-2-3d.tz" known_hash = "6a5de26466028c673d253ca014c75c719467ec6c28d7178baf9287b44ad15191" -downloaded_file_path = pooch.retrieve( - url=URL, known_hash=known_hash, progressbar=True -) +downloaded_file_path = pooch.retrieve(url=URL, known_hash=known_hash, progressbar=True) result = ndk.load_result_from_disk(downloaded_file_path) # %% @@ -32,21 +32,28 @@ # pip install "napari[all]" # ``` # -# or by following the `napari` [installation instructions](https://napari.org/stable/tutorials/fundamentals/installation.html). +# or by following the `napari` installation instructions +# [link](https://napari.org/stable/tutorials/fundamentals/installation.html). try: - import napari + import napari # noqa: F401 + + assert isinstance(result, ndk.results.SteadyStateResult3D) result.render_steady_state_amplitudes_3d() except ImportError: - print("napari has not been installed. Please install it with: pip install napari[all]") + print( + "napari has not been installed. Please install it with: pip install napari[all]" + ) # %% # If you have napari installed you should see an output like the following: # # ``` -# Opening the napari viewer. The window might not show up on top of your notebook; look through your open applications if it does not.""" +# Opening the napari viewer. The window might not show up on top of your notebook; +# look through your open applications if it does not. # ``` # %% -# If you have napari installed you should have been able to see an image like the following: -# ![3d-visualization](../../images/3d_visualization.gif) \ No newline at end of file +# If you have napari installed you should have been able to see an image like +# the following: +# ![3d-visualization](../../images/3d_visualization.gif) diff --git a/docs/examples/plot_adding_custom_source.py b/docs/examples/plot_adding_custom_source.py index 22615213..39a7bd1e 100644 --- a/docs/examples/plot_adding_custom_source.py +++ b/docs/examples/plot_adding_custom_source.py @@ -5,7 +5,8 @@ ==================================== !!! note - NDK and its examples are under constant development, more information and content will be added to this example soon! + NDK and its examples are under constant development, more information + and content will be added to this example soon! This example demonstrates how to add a source to the simulation. """ @@ -27,6 +28,7 @@ # emitting. Defaults to 0.0. import numpy as np + import neurotechdevkit as ndk source = ndk.sources.FocusedSource2D( @@ -37,9 +39,10 @@ num_points=1000, ) -scenario = ndk.make('scenario-0-v0') +scenario = ndk.make("scenario-0-v0") scenario.add_source(source) result = scenario.simulate_steady_state() +assert isinstance(result, ndk.results.SteadyStateResult2D) result.render_steady_state_amplitudes() # %% diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py new file mode 100644 index 00000000..a77bad68 --- /dev/null +++ b/docs/examples/plot_customized_center_frequency.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*- +""" +Custom center frequency +==================================== + +This example demonstrates how to use a customized center frequency +using ndk +""" +# %% +import neurotechdevkit as ndk + +CENTER_FREQUENCY = 6e5 + +scenario = ndk.make("scenario-0-v0") + +# using default material layers +scenario.material_layers = ["water", "cortical_bone", "brain", "tumor"] + +# Customizing material properties +scenario.material_properties = { + "tumor": ndk.materials.Material( + vp=1850.0, rho=1250.0, alpha=0.8, render_color="#94332F" + ), +} + +result = scenario.simulate_steady_state(center_frequency=CENTER_FREQUENCY) +assert isinstance(result, ndk.results.SteadyStateResult2D) +result.render_steady_state_amplitudes(show_material_outlines=False) + +# %% diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py index 1e75d31b..6cc4e0d3 100644 --- a/docs/examples/plot_full_scenario.py +++ b/docs/examples/plot_full_scenario.py @@ -4,7 +4,8 @@ ==================================== !!! note - NDK and its examples are under constant development, more information and content will be added to this example soon! + NDK and its examples are under constant development, more information + and content will be added to this example soon! The following code is a simplified implementation of NDK's Scenario 1. """ @@ -13,10 +14,14 @@ import numpy as np import stride -from mosaic.types import Struct - from neurotechdevkit import sources -from neurotechdevkit.scenarios import Scenario2D, Target, make_grid, add_material_fields_to_problem, materials +from neurotechdevkit.results import SteadyStateResult2D +from neurotechdevkit.scenarios import ( + Scenario2D, + Target, + add_material_fields_to_problem, + make_grid, +) class FullScenario(Scenario2D): @@ -28,6 +33,7 @@ class FullScenario(Scenario2D): doi: 10.1121/10.0013426 https://asa.scitation.org/doi/pdf/10.1121/10.0013426 """ + _SCENARIO_ID = "the_id_for_this_scenario" """ @@ -41,6 +47,17 @@ class FullScenario(Scenario2D): "target_1": Target("target_1", np.array([0.064, 0.0]), 0.004, ""), } + """ + The order of returned materials defines the layering of the scenario. + """ + material_layers = [ + "water", + "skin", + "cortical_bone", + "trabecular_bone", + "brain", + ] + def __init__(self, complexity="fast"): """ Instantiate a new scenario. @@ -53,27 +70,6 @@ def __init__(self, complexity="fast"): complexity=complexity, ) - @property - def _material_layers(self) -> list[tuple[str, Struct]]: - """" - This function defines the material layers for the scenario. - The order of returned materials defines the layering of the scenario. - - Each material has the following properties: - - vp: the speed of sound (in m/s). - - rho: the mass density (in kg/m³). - - alpha: the absorption (in dB/cm). - - render_color: the color used when rendering this material in the scenario layout - plot. - """ - return [ - ("water", materials.water), - ("skin", materials.skin), - ("cortical bone", materials.cortical_bone), - ("trabecular bone", materials.trabecular_bone), - ("brain", materials.brain), - ] - @property def _material_outline_upsample_factor(self) -> int: """ @@ -83,18 +79,17 @@ def _material_outline_upsample_factor(self) -> int: """ return 8 - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency) -> stride.Problem: """The problem definition for the scenario.""" extent = np.array([0.12, 0.07]) # m # scenario constants speed_water = 1500 # m/s - c_freq = 500e3 # hz # desired resolution for complexity=fast ppw = 6 # compute resolution - dx = speed_water / c_freq / ppw # m + dx = speed_water / center_frequency / ppw # m grid = make_grid(extent=extent, dx=dx) problem = stride.Problem( @@ -102,11 +97,11 @@ def _compile_problem(self) -> stride.Problem: ) problem = add_material_fields_to_problem( problem=problem, - materials=self.materials, + materials=self.get_materials(center_frequency), layer_ids=self.layer_ids, masks={ name: _create_scenario_1_mask(name, problem.grid) - for name in self.materials.keys() + for name in self.material_layers }, ) return problem @@ -146,11 +141,11 @@ def _create_scenario_1_mask(material, grid): elif material == "skin": _fill_mask(mask, start=interfaces[0], end=interfaces[1], dx=dx) - elif material == "cortical bone": + elif material == "cortical_bone": _fill_mask(mask, start=interfaces[1], end=interfaces[2], dx=dx) _fill_mask(mask, start=interfaces[3], end=interfaces[4], dx=dx) - elif material == "trabecular bone": + elif material == "trabecular_bone": _fill_mask(mask, start=interfaces[2], end=interfaces[3], dx=dx) elif material == "brain": @@ -178,4 +173,5 @@ def _fill_mask(mask, start, end, dx): scenario = FullScenario() result = scenario.simulate_steady_state() +assert isinstance(result, SteadyStateResult2D) result.render_steady_state_amplitudes(show_material_outlines=False) diff --git a/docs/examples/plot_metrics.py b/docs/examples/plot_metrics.py index 7682a727..7ccc4481 100644 --- a/docs/examples/plot_metrics.py +++ b/docs/examples/plot_metrics.py @@ -5,7 +5,8 @@ ==================================== !!! note - NDK and its examples are under constant development, more information and content will be added to this example soon! + NDK and its examples are under constant development, more information + and content will be added to this example soon! This example demonstrates how to display the metrics collected from the simulation. """ @@ -13,9 +14,9 @@ # ## Rendering scenario import neurotechdevkit as ndk - -scenario = ndk.make('scenario-0-v0') +scenario = ndk.make("scenario-0-v0") result = scenario.simulate_steady_state() +assert isinstance(result, ndk.results.SteadyStateResult2D) result.render_steady_state_amplitudes(show_material_outlines=False) # %% @@ -23,5 +24,7 @@ for metric, metric_value in result.metrics.items(): print(f"{metric}:") print(f"\t {metric_value['description']}") - print(f"\t Unit: {metric_value['unit-of-measurement']} Value: {metric_value['value']}") + print( + f"\t Unit: {metric_value['unit-of-measurement']} Value: {metric_value['value']}" + ) print() diff --git a/docs/examples/plot_multiple_sources.py b/docs/examples/plot_multiple_sources.py index 3bd074af..f4a0def1 100644 --- a/docs/examples/plot_multiple_sources.py +++ b/docs/examples/plot_multiple_sources.py @@ -4,7 +4,8 @@ ==================================== !!! note - NDK and its examples are under constant development, more information and content will be added to this example soon! + NDK and its examples are under constant development, more information + and content will be added to this example soon! Adding multiple sources for transcranial ultrasound stimulation enables greater precision and control in targeting specific areas of the brain. @@ -15,31 +16,32 @@ """ # %% import numpy as np + import neurotechdevkit as ndk -scenario = ndk.make('scenario-0-v0') +scenario = ndk.make("scenario-0-v0") s1 = ndk.sources.FocusedSource2D( - position=np.array([0.01, 0.0]), - direction=np.array([0.92, 0.25]), - aperture=0.01, - focal_length=0.022, - num_points=1000, - ) + position=np.array([0.01, 0.0]), + direction=np.array([0.92, 0.25]), + aperture=0.01, + focal_length=0.022, + num_points=1000, +) s2 = ndk.sources.FocusedSource2D( - position=np.array([0.04, -0.002]), - direction=np.array([-0.85, 0.35]), - aperture=0.01, - focal_length=0.011, - num_points=1000, - delay=5.1e-6, - ) + position=np.array([0.04, -0.002]), + direction=np.array([-0.85, 0.35]), + aperture=0.01, + focal_length=0.011, + num_points=1000, + delay=5.1e-6, +) scenario.add_source(s1) scenario.add_source(s2) result = scenario.simulate_pulse() - +assert isinstance(result, ndk.results.PulsedResult2D) result.render_pulsed_simulation_animation() # %% diff --git a/docs/examples/plot_phased_array_source.py b/docs/examples/plot_phased_array_source.py index acc9ce8b..ac44d3eb 100644 --- a/docs/examples/plot_phased_array_source.py +++ b/docs/examples/plot_phased_array_source.py @@ -3,26 +3,55 @@ Phased array source ==================================== -Phased arrays consist of multiple individual transducer elements arranged in a specific pattern, such as a linear or circle. Each element can be controlled independently, allowing for precise manipulation of the generated ultrasound waves. By adjusting the timing or phase of the signals sent to these elements, the ultrasound waves can be focused, steered, and shaped without mechanically moving the transducer. - -They are becoming increasingly popular in the field of transcranial ultrasound stimulation as they offer several advantages over traditional ultrasound transducers -- among others are precise targeting and better control over steering and shaping: - -* **Precise targeting:** as individual transducer elements can be individually controlled, it allows for the generation of a focused ultrasound beam with high spatial accuracy. This enables the stimulation of specific brain regions without affecting the surrounding healthy tissue and minimizes the risk of potential side effects. - -* **Steering and shaping:** The phased array technology allows the ultrasound beam to be electronically steered and shaped in real-time without mechanically moving the transducers. This enables targeting of different brain regions or adjusting the stimulation pattern as needed during a session, making the procedure more versatile and adaptable. - -These features allow the stimulation to be tailored to suit to more specific requirements. With ongoing research and development, they have the potential to revolutionize the field of brain stimulation and offer new treatment options for a range of neurological and psychiatric disorders. - -This notebook will show how to define a phased array within NDK and experiment with some of the available features. For more details, checkout the [NDK documentation](https://agencyenterprise.github.io/neurotechdevkit/). +Phased arrays consist of multiple individual transducer elements arranged in a +specific pattern, such as a linear or circle. Each element can be controlled +independently, allowing for precise manipulation of the generated ultrasound +waves. By adjusting the timing or phase of the signals sent to these elements, +the ultrasound waves can be focused, steered, and shaped without mechanically +moving the transducer. + +They are becoming increasingly popular in the field of transcranial ultrasound +stimulation as they offer several advantages over traditional ultrasound +transducers -- among others are precise targeting and better control over steering +and shaping: + +* **Precise targeting:** as individual transducer elements can be individually +controlled, it allows for the generation of a focused ultrasound beam with high +spatial accuracy. This enables the stimulation of specific brain regions without +affecting the surrounding healthy tissue and minimizes the risk of potential side +effects. + +* **Steering and shaping:** The phased array technology allows the ultrasound beam +to be electronically steered and shaped in real-time without mechanically moving the +transducers. This enables targeting of different brain regions or adjusting the +stimulation pattern as needed during a session, making the procedure more versatile +and adaptable. + +These features allow the stimulation to be tailored to suit to more specific +requirements. With ongoing research and development, they have the potential to +revolutionize the field of brain stimulation and offer new treatment options for +a range of neurological and psychiatric disorders. + +This notebook will show how to define a phased array within NDK and experiment with +some of the available features. For more details, checkout the +[NDK documentation](https://agencyenterprise.github.io/neurotechdevkit/). ## Phased Arrays Components -* **Elements:** A phased array consists of multiple transducer elements that can be individually controlled in terms of phase. -* **Pitch:** The pitch of a phased array refers to the distance between adjacent transducer elements, and affects the resolution of the array. -* **Spacing:** The spacing between transducer elements affects the width and shape of the ultrasound beam that is produced. -* **Aperture:** The aperture of a phased array refers to the total size of the array, and affects the power and focus of the ultrasound beam. -* **Element width:** The element width of a phased array refers to the size of each individual transducer element. Wider elements generally result in a more powerful and focused beam, while narrower elements can increase the resolution of the array. -* **Tilt angle**: The tilt angle of a phased array refers to the angle at which the ultrasound beam is directed relative to the normal of plane of the array. Tilt angle can be adjusted by controlling the phase (delays) of the individual transducer elements. +* **Elements:** A phased array consists of multiple transducer elements that can be +individually controlled in terms of phase. +* **Pitch:** The pitch of a phased array refers to the distance between adjacent +transducer elements, and affects the resolution of the array. +* **Spacing:** The spacing between transducer elements affects the width and shape +of the ultrasound beam that is produced. +* **Aperture:** The aperture of a phased array refers to the total size of the array, +and affects the power and focus of the ultrasound beam. +* **Element width:** The element width of a phased array refers to the size of each +individual transducer element. Wider elements generally result in a more powerful +and focused beam, while narrower elements can increase the resolution of the array. +* **Tilt angle**: The tilt angle of a phased array refers to the angle at which the +ultrasound beam is directed relative to the normal of plane of the array. Tilt angle +can be adjusted by controlling the phase (delays) of the individual transducer elements. ![phased-array-schematic](../../images/phased-array-schematic.png) @@ -31,10 +60,13 @@ ![phased-array-axis-definition](../../images/phased-array-axis-definition.png) -The image above shows the axis definition and the look of a phased array in the real world. Image source [link](https://www.biomecardio.com/MUST/functions/html/txdelay_doc.html). +The image above shows the axis definition and the look of a phased array in the real +world. Image source +[link](https://www.biomecardio.com/MUST/functions/html/txdelay_doc.html). -Below we will go through some examples of how to use the NDK to API to simulate Phased Arrays in the most common settings. Examples we wil cover: +Below we will go through some examples of how to use the NDK to API to simulate +Phased Arrays in the most common settings. Examples we will cover: * Setting up a tilted wavefront * Focusing phased arrays @@ -58,10 +90,11 @@ - num_elements `(int)`: the number of elements of the phased array. - pitch `(float)`: the distance (in meters) between the centers of neighboring elements in the phased array. -- element_width `(float)`: the width (in meters) of each individual element of the array. -- tilt_angle `(float)`: the desired tilt angle (in degrees) of the wavefront. The angle is - measured between the direction the wavefront travels and the normal to the - surface of the transducer, with positive angles resulting in a +- element_width `(float)`: the width (in meters) of each individual element of + the array. +- tilt_angle `(float)`: the desired tilt angle (in degrees) of the wavefront. + The angle is measured between the direction the wavefront travels and the + normal to the surface of the transducer, with positive angles resulting in a counter-clockwise tilt away from the normal. - focal_length `(float)`: the distance (in meters) from `position` to the focal point. @@ -81,75 +114,80 @@ # argument different than zero. Positive angles result in counter-clockwise rotations. import numpy as np -import neurotechdevkit as ndk +import neurotechdevkit as ndk # define the source source = ndk.sources.PhasedArraySource2D( - position=np.array([0.0, 0.0]), - direction=np.array([1., 0.]), - num_elements=20, - pitch=1.5e-3, - element_width=1.2e-3, - tilt_angle=30, - num_points=1000, + position=np.array([0.0, 0.0]), + direction=np.array([1.0, 0.0]), + num_elements=20, + pitch=1.5e-3, + element_width=1.2e-3, + tilt_angle=30, + num_points=1000, ) -scenario = ndk.make('scenario-1-2d-v0') +scenario = ndk.make("scenario-1-2d-v0") scenario.add_source(source) result = scenario.simulate_steady_state() +assert isinstance(result, ndk.results.SteadyStateResult2D) result.render_steady_state_amplitudes() # %% # ## Focusing -# Similarly, a phased array can be focused by applying delays in a specific way. Such delays are automatically computed when `focal_length` is defined. Let's explore how the API looks like: +# Similarly, a phased array can be focused by applying delays in a specific way. Such +# delays are automatically computed when `focal_length` is defined. Let's explore how +# the API looks like: -scenario = ndk.make('scenario-1-2d-v0') +scenario = ndk.make("scenario-1-2d-v0") phased_array = ndk.sources.PhasedArraySource2D( - position=np.array([0.0, 0.0]), - direction=np.array([1.0, 0.0]), - num_elements=20, - pitch=1.5e-3, - element_width=1.2e-3, - focal_length=0.03, - num_points=1000, + position=np.array([0.0, 0.0]), + direction=np.array([1.0, 0.0]), + num_elements=20, + pitch=1.5e-3, + element_width=1.2e-3, + focal_length=0.03, + num_points=1000, ) scenario.add_source(phased_array) -print(f'Focal point is: {phased_array.focal_point}') +print(f"Focal point is: {phased_array.focal_point}") # %% # `focal_point` shows the coordinates (in meters) where the array focuses. result = scenario.simulate_steady_state() +assert isinstance(result, ndk.results.SteadyStateResult2D) result.render_steady_state_amplitudes() # %% # We can see that the arrays focuses at a distance equal to `focal_length`. - # %% # ## Custom Wavefront Shapes -# In the examples shown so far we specified high level arguments and the required delays were -# automatically computed. The NDK API offers also the possibility to specify delays directly -# to create arbitrary wavefront shapes. +# In the examples shown so far we specified high level arguments and the required +# delays were automatically computed. The NDK API offers also the possibility to +# specify delays directly to create arbitrary wavefront shapes. # -# In the example below we apply monotonically increasing delays to mimic a counter-clockwise angle. +# In the example below we apply monotonically increasing delays to mimic a +# counter-clockwise angle. -scenario = ndk.make('scenario-1-2d-v0') +scenario = ndk.make("scenario-1-2d-v0") phased_array = ndk.sources.PhasedArraySource2D( - position=np.array([0.0, 0.0]), - direction=np.array([1., 0.]), - num_elements=20, - pitch=1.5e-3, - element_width=1.2e-3, - element_delays=np.linspace(0, 1e-5, 20), - num_points=1000, + position=np.array([0.0, 0.0]), + direction=np.array([1.0, 0.0]), + num_elements=20, + pitch=1.5e-3, + element_width=1.2e-3, + element_delays=np.linspace(0, 1e-5, 20), + num_points=1000, ) scenario.add_source(phased_array) result = scenario.simulate_steady_state() +assert isinstance(result, ndk.results.SteadyStateResult2D) result.render_steady_state_amplitudes() # %% @@ -166,21 +204,23 @@ # # The rest of the API is identical for both 2D and 3D scenarios. -scenario = ndk.make('scenario-1-3d-v0') +scenario = ndk.make("scenario-1-3d-v0") phased_3d = ndk.sources.PhasedArraySource3D( - position=np.array([0.,0.,0.]), - direction=np.array([1.,0.,0.]), - center_line=np.array([0.,0.,1.]), + position=np.array([0.0, 0.0, 0.0]), + direction=np.array([1.0, 0.0, 0.0]), + center_line=np.array([0.0, 0.0, 1.0]), num_points=20_000, num_elements=16, pitch=1.5e-3, element_width=1.2e-3, tilt_angle=30, - height=5.e-3) + height=5.0e-3, +) scenario.add_source(phased_3d) results = scenario.simulate_steady_state() +assert isinstance(results, ndk.results.SteadyStateResult3D) results.render_steady_state_amplitudes(slice_axis=1, slice_position=0.0) # %% diff --git a/docs/examples/plot_pulsed_simulation.py b/docs/examples/plot_pulsed_simulation.py index 067d2f59..dee1c005 100644 --- a/docs/examples/plot_pulsed_simulation.py +++ b/docs/examples/plot_pulsed_simulation.py @@ -11,13 +11,16 @@ scenario = ndk.make("scenario-0-v0") result = scenario.simulate_pulse() +assert isinstance(result, ndk.results.PulsedResult2D) result.render_pulsed_simulation_animation() # %% # # Generating a video file # ======================= -# You can also generate a video file of the simulation (which requires [ffmpeg](https://ffmpeg.org/download.html) installed). +# You can also generate a video file of the simulation (which requires +# [ffmpeg](https://ffmpeg.org/download.html) installed). # -# To create and save the video as `animation.mp4` in the current folder, all you need is the execute following command: +# To create and save the video as `animation.mp4` in the current folder, all you +# need is the execute following command: result.create_video_file("animation.mp4", fps=25, overwrite=True) diff --git a/docs/examples/plot_scenarios.py b/docs/examples/plot_scenarios.py index 1711704c..fcdbadb5 100644 --- a/docs/examples/plot_scenarios.py +++ b/docs/examples/plot_scenarios.py @@ -11,22 +11,23 @@ def plot_scenario(scenario_id): - print(f"Simulating scenario: {scenario_id}") - scenario = ndk.make(scenario_id) - result = scenario.simulate_steady_state() - result.render_steady_state_amplitudes(show_material_outlines=False) + print(f"Simulating scenario: {scenario_id}") + scenario = ndk.make(scenario_id) + result = scenario.simulate_steady_state() + result.render_steady_state_amplitudes(show_material_outlines=False) + # %% # Simulating scenario: scenario-0-v0 # =================================== -plot_scenario('scenario-0-v0') +plot_scenario("scenario-0-v0") # %% # Simulating scenario: scenario-1-2d-v0 # =================================== -plot_scenario('scenario-1-2d-v0') +plot_scenario("scenario-1-2d-v0") # %% # Simulating scenario: scenario-2-2d-v0 # =================================== -plot_scenario('scenario-2-2d-v0') +plot_scenario("scenario-2-2d-v0") diff --git a/docs/examples/plot_store_results.py b/docs/examples/plot_store_results.py index 51fc0403..bafdec6a 100644 --- a/docs/examples/plot_store_results.py +++ b/docs/examples/plot_store_results.py @@ -4,26 +4,29 @@ ==================================== !!! note - NDK and its examples are under constant development, more information and content will be added to this example soon! + NDK and its examples are under constant development, more information and + content will be added to this example soon! This example demonstrates how a simulation can be executed and stored in one -computer and the results can be loaded and rendered later in the same computer or another one. +computer and the results can be loaded and rendered later in the same computer +or another one. """ # %% # Execute the following code in a computer with ndk installed import neurotechdevkit as ndk -scenario = ndk.make('scenario-0-v0') +scenario = ndk.make("scenario-0-v0") result = scenario.simulate_steady_state() -result.save_to_disk('scenario-0-v0-results.tar.gz') +result.save_to_disk("scenario-0-v0-results.tar.gz") # %% # The output file from the previous step could be copied to # another computer with ndk installed or stored for future use. The results can # be loaded running the following code: -import neurotechdevkit as ndk +import neurotechdevkit as ndk # noqa: E402 -result = ndk.load_result_from_disk('scenario-0-v0-results.tar.gz') -result.render_steady_state_amplitudes(show_material_outlines=False) +result2 = ndk.load_result_from_disk("scenario-0-v0-results.tar.gz") +assert isinstance(result2, ndk.results.SteadyStateResult2D) +result2.render_steady_state_amplitudes(show_material_outlines=False) # %% diff --git a/src/neurotechdevkit/__init__.py b/src/neurotechdevkit/__init__.py index 3cefd11c..14b97161 100644 --- a/src/neurotechdevkit/__init__.py +++ b/src/neurotechdevkit/__init__.py @@ -3,12 +3,13 @@ import os -from . import scenarios, sources +from . import materials, scenarios, sources from .results import load_result_from_disk __all__ = [ "results", "scenarios", + "materials", "sources", "make", "ScenarioNotFoundError", diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py new file mode 100644 index 00000000..805e74ec --- /dev/null +++ b/src/neurotechdevkit/materials.py @@ -0,0 +1,157 @@ +"""Materials for the neurotechdevkit scenarios.""" +from dataclasses import dataclass + +from mosaic.types import Struct + +NEPER_TO_DECIBEL = 8.6860000037 + + +@dataclass +class _BaseMaterial: + """A material with properties used in the neurotechdevkit scenarios.""" + + vp: float + rho: float + render_color: str + + +@dataclass +class Material(_BaseMaterial): + """A NDK Material with an attenuation coefficient.""" + + alpha: float + + def to_struct(self) -> Struct: + """Return a Struct representation of the material. + + Returns: + Struct: a Struct representation of the material. + """ + struct = Struct() + struct.vp = self.vp + struct.rho = self.rho + struct.alpha = self.alpha + struct.render_color = self.render_color + return struct + + +# The values below consider a center frequency of 500kHz +DEFAULT_MATERIALS = { + "water": _BaseMaterial(vp=1500.0, rho=1000.0, render_color="#2E86AB"), + "skin": _BaseMaterial(vp=1610.0, rho=1090.0, render_color="#FA8B53"), + "cortical_bone": _BaseMaterial(vp=2800.0, rho=1850.0, render_color="#FAF0CA"), + "trabecular_bone": _BaseMaterial(vp=2300.0, rho=1700.0, render_color="#EBD378"), + "brain": _BaseMaterial(vp=1560.0, rho=1040.0, render_color="#DB504A"), + # these numbers are completely made up + # TODO: research reasonable values + "tumor": _BaseMaterial(vp=1650.0, rho=1150.0, render_color="#94332F"), +} + + +@dataclass +class AttenuationConstant: + """The parameters of the attenuation constant.""" + + a0: float # α0 [Np/m/MHz] + b: float + + def calculate_absorption(self, frequency: float) -> float: + """ + Calculate the absorption coefficient for a given center frequency. + + The absorption coefficient is calculated using the following formula: + α(f) = α0 * f^b + + We convert the absorption coefficient from Np/m/Hz to dB/m/Hz + multiplying it by 8.6860000037. + + Args: + frequency (float): in Hz + + Returns: + float: attenuation in dB/cm/MHz + """ + attenuation = self.a0 * (frequency / 1e6) ** self.b # Np/m + db_attenuation = attenuation * NEPER_TO_DECIBEL # dB/m + return db_attenuation / 100 # convert from dB/m to dB/cm + + +ATTENUATION_CONSTANTS = { + "trabecular_bone": AttenuationConstant(a0=47, b=1.2), + "cortical_bone": AttenuationConstant(a0=54.553, b=1), + "skin": AttenuationConstant(a0=21.158, b=1), + "water": AttenuationConstant(a0=0.025328436, b=1), + "brain": AttenuationConstant(a0=6.8032, b=1.3), + # these numbers are completely made up + # TODO: research reasonable values + "tumor": AttenuationConstant(a0=8, b=1), +} + + +def _calculate_absorption(material_name: str, center_frequency: float) -> float: + """ + Calculate the absorption coefficient for a given center frequency. + + Args: + material_name (str): the name of the material. E.g. "water" + center_frequency (float): the center frequency of the transducer in Hz + + Raises: + ValueError: raised if the material is not supported. + + Returns: + float: the absorption coefficient in dB/cm/Hz + """ + try: + return ATTENUATION_CONSTANTS[material_name].calculate_absorption( + center_frequency + ) + except KeyError: + raise ValueError(f"Undefined material: {material_name}") + + +def get_render_color(material_name: str) -> str: + """ + Get the render color for a material. + + Args: + material_name (str): the name of the material. E.g. "water" + + Raises: + ValueError: raised if the material is not supported. + + Returns: + str: the render color of the material. + """ + if material_name not in DEFAULT_MATERIALS: + raise ValueError(f"Undefined material: {material_name}") + + return DEFAULT_MATERIALS[material_name].render_color + + +def get_material(material_name: str, center_frequency: float = 5.0e5) -> Material: + """ + Get a material with properties used in the neurotechdevkit scenarios. + + Args: + material_name (str): the name of the material. E.g. "water" + center_frequency (float, optional): The center frequency of the + transducer. Defaults to 5.0e5. + + Raises: + ValueError: raised if the material is not supported. + + Returns: + Material: a material with properties used in the neurotechdevkit. + """ + if material_name not in DEFAULT_MATERIALS: + raise ValueError(f"Undefined material: {material_name}") + + base_material = DEFAULT_MATERIALS[material_name] + + return Material( + vp=base_material.vp, + rho=base_material.rho, + alpha=_calculate_absorption(material_name, center_frequency), + render_color=base_material.render_color, + ) diff --git a/src/neurotechdevkit/rendering/napari.py b/src/neurotechdevkit/rendering/napari.py index 2d049506..8c7af3b4 100644 --- a/src/neurotechdevkit/rendering/napari.py +++ b/src/neurotechdevkit/rendering/napari.py @@ -98,7 +98,6 @@ def render_amplitudes_3d_with_napari(result: "results.SteadyStateResult3D") -> N Raises: ImportError: If napari is not found. """ - pass _create_napari_3d(scenario=result.scenario, amplitudes=result.get_steady_state()) diff --git a/src/neurotechdevkit/results/_metrics.py b/src/neurotechdevkit/results/_metrics.py index 63ec8813..ece788d7 100644 --- a/src/neurotechdevkit/results/_metrics.py +++ b/src/neurotechdevkit/results/_metrics.py @@ -113,7 +113,7 @@ def calculate_all_metrics( " intended use)." ), } - for layer in result.scenario.ordered_layers + for layer in result.scenario.material_layers }, } diff --git a/src/neurotechdevkit/results/_results.py b/src/neurotechdevkit/results/_results.py index 0f091e6b..b0a5b574 100644 --- a/src/neurotechdevkit/results/_results.py +++ b/src/neurotechdevkit/results/_results.py @@ -1297,6 +1297,7 @@ def load_result_from_disk(filepath: str | pathlib.Path) -> Result: wavefield=save_data.get("wavefield"), traces=None, ) + scenario._problem = scenario._compile_problem(save_data["center_frequency"]) if save_data.get("steady_state") is not None: fields_kwargs.update(steady_state=save_data["steady_state"]) diff --git a/src/neurotechdevkit/scenarios/__init__.py b/src/neurotechdevkit/scenarios/__init__.py index dc4d0181..d10110e5 100644 --- a/src/neurotechdevkit/scenarios/__init__.py +++ b/src/neurotechdevkit/scenarios/__init__.py @@ -1,5 +1,5 @@ """Scenarios module.""" -from . import materials +from .. import materials from ._base import Scenario, Scenario2D, Scenario3D, Target from ._scenario_0 import Scenario0 from ._scenario_1 import Scenario1_2D, Scenario1_3D diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py index 2e0c28e0..3949482c 100644 --- a/src/neurotechdevkit/scenarios/_base.py +++ b/src/neurotechdevkit/scenarios/_base.py @@ -16,6 +16,7 @@ from stride.problem import StructuredData from .. import rendering, results +from ..materials import Material, get_material, get_render_color from ..sources import Source from ._resources import budget_time_and_memory_resources from ._shots import create_shot @@ -60,6 +61,12 @@ class Scenario(abc.ABC): _SCENARIO_ID: str _TARGET_OPTIONS: dict[str, Target] + # The list of material layers in the scenario. + material_layers: list[str] + + # The customization to the material layers. + material_properties: dict[str, Material] = {} + def __init__( self, origin: npt.NDArray[np.float_], @@ -71,7 +78,6 @@ def __init__( raise ValueError("the only complexity currently supported is 'fast'") self._origin = origin - self._problem = self._compile_problem() self._sources: FrozenList[Source] = FrozenList() self._target_id: str @@ -214,9 +220,8 @@ def target_radius(self) -> float: """The radius of the target region (in meters).""" return self.target.radius - @property - def materials(self) -> Mapping[str, Struct]: - """A map between material name and material properties. + def get_materials(self, center_frequency=float) -> Mapping[str, Struct]: + """Return a map between material name and material properties. - vp: the speed of sound (in m/s). - rho: the mass density (in kg/m³). @@ -224,22 +229,36 @@ def materials(self) -> Mapping[str, Struct]: - render_color: the color used when rendering this material in the scenario layout plot. """ - return {name: material for name, material in self._material_layers} + materials = {} + for layer in self.material_layers: + if layer not in self.material_properties: + material_properties = get_material(layer, center_frequency) + else: + material_properties = self.material_properties[layer] + materials[layer] = material_properties.to_struct() + return materials @property - def layer_ids(self) -> Mapping[str, int]: - """A map between material names and their layer id.""" - return {name: n for n, (name, _) in enumerate(self._material_layers)} + def material_colors(self) -> dict[str, str]: + """ + A map between material name and material render color. - @property - def ordered_layers(self) -> list[str]: - """A list of material names in order of their layer id.""" - return [name for name, _ in self._material_layers] + Returns: + dict[str, str]: keys are material names and values are the hex color + """ + material_colors = {} + for material in self.material_layers: + if material in self.material_properties: + color = self.material_properties[material].render_color + else: + color = get_render_color(material_name=material) + material_colors[material] = color + return material_colors @property - @abc.abstractmethod - def _material_layers(self) -> list[tuple[str, Struct]]: - pass + def layer_ids(self) -> Mapping[str, int]: + """A map between material names and their layer id.""" + return {name: n for n, name in enumerate(self.material_layers)} @property @abc.abstractmethod @@ -320,7 +339,7 @@ def get_field_data(self, field: str) -> npt.NDArray[np.float_]: return self.problem.medium.fields[field].data @abc.abstractmethod - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: pass def reset(self) -> None: @@ -372,10 +391,6 @@ def simulate_steady_state( found by taking the Fourier transform of the last `n_cycles_steady_state` cycles of data and taking the amplitude of the component at the `center_frequency`. - !!! note - The only supported frequency currently supported is 500kHz. Any other - value will raise a NotImplementedError. - !!! warning A poor choice of arguments to this function can lead to a failed simulation. Make sure you understand the impact of supplying parameter @@ -383,8 +398,7 @@ def simulate_steady_state( Args: center_frequency: the center frequency (in hertz) to use for the - continuous-wave source output. No other value besides 500kHz (the - default) is currently supported. + continuous-wave source output. points_per_period: the number of points in time to simulate for each cycle of the wave. n_cycles_steady_state: the number of complete cycles to use when calculating @@ -400,23 +414,14 @@ def simulate_steady_state( n_jobs: the number of threads to be used for the computation. Use None to leverage Devito automatic tuning. - Raises: - NotImplementedError: if a `center_frequency` other than 500kHz is provided. - Returns: An object containing the result of the steady-state simulation. """ - if center_frequency != 5.0e5: - raise NotImplementedError( - "500kHz is the only currently supported center frequency. Support for" - " other frequencies will be implemented once material properties as a" - " function of frequency has been implemented." - ) - + self._problem = self._compile_problem(center_frequency) problem = self.problem sim_time = select_simulation_time_for_steady_state( grid=problem.grid, - materials=self.materials, + materials=self.get_materials(center_frequency), freq_hz=center_frequency, time_to_steady_state=time_to_steady_state, n_cycles_steady_state=n_cycles_steady_state, @@ -474,10 +479,6 @@ def simulate_pulse( In this simulation, the sources will emit a pulse containing a few cycles of oscillation and then let the pulse propagate out to all edges of the scenario. - !!! note - The only supported frequency currently supported is 500kHz. Any other - value will raise a NotImplementedError. - !!! warning A poor choice of arguments to this function can lead to a failed simulation. Make sure you understand the impact of supplying parameter @@ -485,8 +486,7 @@ def simulate_pulse( Args: center_frequency: the center frequency (in hertz) to use for the - continuous-wave source output. No other value besides - 500kHz (the default) is currently supported. + continuous-wave source output. points_per_period: the number of points in time to simulate for each cycle of the wave. simulation_time: the amount of time (in seconds) the simulation should run. @@ -499,9 +499,6 @@ def simulate_pulse( n_jobs: the number of threads to be used for the computation. Use None to leverage Devito automatic tuning. - Raises: - NotImplementedError: if a `center_frequency` other than 500kHz is provided. - Returns: An object containing the result of the 2D pulsed simulation. """ @@ -530,18 +527,13 @@ def _simulate_pulse( In this simulation, the sources will emit a pulse containing a few cycles of oscillation and then let the pulse propagate out to all edges of the scenario. - !!! note - The only supported frequency currently supported is 500kHz. Any other - value will raise a NotImplementedError. - Warning: A poor choice of arguments to this function can lead to a failed simulation. Make sure you understand the impact of supplying parameter values other than the default if you chose to do so. Args: center_frequency: the center frequency (in hertz) to use for the - continuous-wave source output. No other value besides - 500kHz (the default) is currently supported. + continuous-wave source output. points_per_period: the number of points in time to simulate for each cycle of the wave. simulation_time: the amount of time (in seconds) the simulation should run. @@ -560,24 +552,15 @@ def _simulate_pulse( which the slice of the 3D field should be made. Only valid if `slice_axis` is not None. - Raises: - NotImplementedError: if a `center_frequency` other than 500kHz is provided. - Returns: An object containing the result of the pulsed simulation. """ - if center_frequency != 5.0e5: - raise NotImplementedError( - "500kHz is the only currently supported center frequency. Support for" - " other frequencies will be implemented once material properties as a" - " function of frequency has been implemented." - ) - + self._problem = self._compile_problem(center_frequency) problem = self.problem if simulation_time is None: simulation_time = select_simulation_time_for_pulsed( grid=problem.grid, - materials=self.materials, + materials=self.get_materials(center_frequency), delay=find_largest_delay_in_sources(self.sources), ) problem.grid.time = create_time_grid( @@ -904,9 +887,7 @@ def render_layout( show_material_outlines: whether or not to display a thin white outline of the transition between different materials. """ - color_sequence = [ - self.materials[name].render_color for name in self.ordered_layers - ] + color_sequence = list(self.material_colors.values()) field = self.get_field_data("layer").astype(int) fig, ax = rendering.create_layout_fig( self.extent, self.origin, color_sequence, field @@ -939,7 +920,7 @@ def render_layout( fig=fig, ax=ax, color_sequence=color_sequence, - layer_labels=self.ordered_layers, + layer_labels=self.material_layers, show_sources=show_sources, show_target=show_target, extent=self.extent, @@ -1008,10 +989,6 @@ def simulate_pulse( In this simulation, the sources will emit a pulse containing a few cycles of oscillation and then let the pulse propagate out to all edges of the scenario. - !!! note - The only supported frequency currently supported is 500kHz. Any - other value will raise a NotImplementedError. - !!! warning A poor choice of arguments to this function can lead to a failed simulation. Make sure you understand the impact of supplying parameter @@ -1019,8 +996,7 @@ def simulate_pulse( Args: center_frequency: the center frequency (in hertz) to use for the - continuous-wave source output. No other value besides - 500kHz (the default) is currently supported. + continuous-wave source output. points_per_period: the number of points in time to simulate for each cycle of the wave. simulation_time: the amount of time (in seconds) the simulation should run. @@ -1039,9 +1015,6 @@ def simulate_pulse( which the slice of the 3D field should be made. Only valid if `slice_axis` is not None. - Raises: - NotImplementedError: if a `center_frequency` other than 500kHz is provided. - Returns: An object containing the result of the 3D pulsed simulation. """ @@ -1087,9 +1060,7 @@ def render_layout( if slice_position is None: slice_position = self.get_default_slice_position(slice_axis) - color_sequence = [ - self.materials[name].render_color for name in self.ordered_layers - ] + color_sequence = list(self.material_colors.values()) field = self.get_field_data("layer").astype(int) field = slice_field(field, self, slice_axis, slice_position) extent = drop_element(self.extent, slice_axis) @@ -1127,7 +1098,7 @@ def render_layout( fig=fig, ax=ax, color_sequence=color_sequence, - layer_labels=self.ordered_layers, + layer_labels=self.material_layers, show_sources=show_sources, show_target=show_target, extent=extent, diff --git a/src/neurotechdevkit/scenarios/_scenario_0.py b/src/neurotechdevkit/scenarios/_scenario_0.py index 2c631ca4..ad0cfc9e 100644 --- a/src/neurotechdevkit/scenarios/_scenario_0.py +++ b/src/neurotechdevkit/scenarios/_scenario_0.py @@ -1,9 +1,8 @@ import numpy as np import stride -from mosaic.types import Struct +from ..materials import Material from ..sources import FocusedSource2D -from . import materials from ._base import Scenario2D, Target from ._utils import ( add_material_fields_to_problem, @@ -25,6 +24,20 @@ class Scenario0(Scenario2D): description="Represents a simulated tumor.", ), } + material_layers = [ + "water", + "cortical_bone", + "brain", + "tumor", + ] + material_properties = { + "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"), + "cortical_bone": Material( + vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA" + ), + "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"), + "tumor": Material(vp=1650.0, rho=1150.0, alpha=0.8, render_color="#94332F"), + } def __init__(self, complexity="fast"): """Create a new instance of scenario 0.""" @@ -35,15 +48,6 @@ def __init__(self, complexity="fast"): complexity=complexity, ) - @property - def _material_layers(self) -> list[tuple[str, Struct]]: - return [ - ("water", materials.water), - ("skull", materials.cortical_bone), - ("brain", materials.brain), - ("tumor", materials.tumor), - ] - @property def _material_outline_upsample_factor(self) -> int: return 16 @@ -51,22 +55,21 @@ def _material_outline_upsample_factor(self) -> int: def _get_material_masks(self, problem): return { name: _create_scenario_0_mask(name, problem.grid, self._origin) - for name in self.materials.keys() + for name in self.material_layers } - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency=float) -> stride.Problem: extent = np.array([0.05, 0.04]) # m origin = self.origin # m # scenario constants speed_water = 1500 # m/s - c_freq = 500e3 # hz # desired resolution for complexity=fast ppw = 6 # compute resolution - dx = speed_water / c_freq / ppw # m + dx = speed_water / center_frequency / ppw # m grid = make_grid(extent=extent, dx=dx) self._origin = origin @@ -75,7 +78,7 @@ def _compile_problem(self) -> stride.Problem: ) problem = add_material_fields_to_problem( problem=problem, - materials=self.materials, + materials=self.get_materials(center_frequency), layer_ids=self.layer_ids, masks=self._get_material_masks(problem), ) @@ -98,7 +101,7 @@ def _create_scenario_0_mask(material, grid, origin): water_mask = ~outer_skull_mask return water_mask - elif material == "skull": + elif material == "cortical_bone": outer_skull_mask = _create_skull_interface_mask(grid, origin) outer_brain_mask = _create_brain_interface_mask(grid, origin) skull_mask = outer_skull_mask & ~outer_brain_mask diff --git a/src/neurotechdevkit/scenarios/_scenario_1.py b/src/neurotechdevkit/scenarios/_scenario_1.py index c728d797..b1fd15ff 100644 --- a/src/neurotechdevkit/scenarios/_scenario_1.py +++ b/src/neurotechdevkit/scenarios/_scenario_1.py @@ -1,45 +1,19 @@ from __future__ import annotations -from typing import Mapping, Protocol +from typing import Mapping import numpy as np import numpy.typing as npt import stride -from mosaic.types import Struct from .. import rendering, sources -from . import materials -from ._base import Scenario2D, Scenario3D, Target +from ..materials import Material +from ._base import Scenario, Scenario2D, Scenario3D, Target from ._utils import add_material_fields_to_problem, make_grid -class _Scenario1MixinProtocol(Protocol): - """Provide type-hinting for Scenario 1 members used by mixins.""" - - @property - def scenario_id(self) -> str: - ... - - @property - def complexity(self) -> str: - ... - - @property - def materials(self) -> Mapping[str, Struct]: - ... - - @property - def layer_ids(self) -> Mapping[str, int]: - ... - - def _get_material_masks( - self, problem: stride.Problem - ) -> Mapping[str, npt.NDArray[np.bool_]]: - ... - - -class _Scenario1Mixin: - """A mixin providing specific implementation detail for scenario 1. +class Scenario1(Scenario): + """Specific implementation detail for scenario 1. Scenario 1 is based on benchmark 4 of the following paper: @@ -51,40 +25,48 @@ class _Scenario1Mixin: https://asa.scitation.org/doi/pdf/10.1121/10.0013426 """ - @property - def _material_layers(self: _Scenario1MixinProtocol) -> list[tuple[str, Struct]]: - return [ - ("water", materials.water), - ("skin", materials.skin), - ("cortical bone", materials.cortical_bone), - ("trabecular bone", materials.trabecular_bone), - ("brain", materials.brain), - ] + material_layers = [ + "water", + "skin", + "cortical_bone", + "trabecular_bone", + "brain", + ] + material_properties = { + "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"), + "skin": Material(vp=1610.0, rho=1090.0, alpha=0.2, render_color="#FA8B53"), + "cortical_bone": Material( + vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA" + ), + "trabecular_bone": Material( + vp=2300.0, rho=1700.0, alpha=8.0, render_color="#EBD378" + ), + "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"), + } @property def _material_outline_upsample_factor(self) -> int: return 8 def _get_material_masks( - self: _Scenario1MixinProtocol, problem: stride.Problem + self, problem: stride.Problem ) -> Mapping[str, npt.NDArray[np.bool_]]: return { name: _create_scenario_1_mask(name, problem.grid) - for name in self.materials.keys() + for name in self.material_layers } def _compile_scenario_1_problem( - self: _Scenario1MixinProtocol, extent: npt.NDArray[np.float_] + self, extent: npt.NDArray[np.float_], center_frequency: float ) -> stride.Problem: # scenario constants speed_water = 1500 # m/s - c_freq = 500e3 # hz # desired resolution for complexity=fast ppw = 6 # compute resolution - dx = speed_water / c_freq / ppw # m + dx = speed_water / center_frequency / ppw # m grid = make_grid(extent=extent, dx=dx) problem = stride.Problem( @@ -92,14 +74,14 @@ def _compile_scenario_1_problem( ) problem = add_material_fields_to_problem( problem=problem, - materials=self.materials, + materials=self.get_materials(center_frequency), layer_ids=self.layer_ids, masks=self._get_material_masks(problem), ) return problem -class Scenario1_2D(_Scenario1Mixin, Scenario2D): +class Scenario1_2D(Scenario1, Scenario2D): """A 2D implementation of scenario 1. Scenario 1 is based on benchmark 4 of the following paper: @@ -126,9 +108,9 @@ def __init__(self, complexity="fast"): complexity=complexity, ) - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: extent = np.array([0.12, 0.07]) # m - return self._compile_scenario_1_problem(extent) + return self._compile_scenario_1_problem(extent, center_frequency) def get_default_source(self) -> sources.Source: """Return the default source for the scenario.""" @@ -141,7 +123,7 @@ def get_default_source(self) -> sources.Source: ) -class Scenario1_3D(_Scenario1Mixin, Scenario3D): +class Scenario1_3D(Scenario1, Scenario3D): """A 3D implementation of scenario 1. Scenario 1 is based on benchmark 4 of the following paper: @@ -190,15 +172,15 @@ def viewer_config_3d(self) -> rendering.ViewerConfig3D: colormaps={ "water": "blue", "skin": "viridis", - "cortical bone": "magma", - "trabecular bone": "inferno", + "cortical_bone": "magma", + "trabecular_bone": "inferno", "brain": "bop orange", }, opacities={ "water": 0.8, "skin": 0.2, - "cortical bone": 0.2, - "trabecular bone": 0.2, + "cortical_bone": 0.2, + "trabecular_bone": 0.2, "brain": 0.4, }, ) @@ -212,9 +194,9 @@ def get_default_slice_position(self, axis: int) -> float: default_positions = np.array([0.064, 0.0, 0.0]) return default_positions[axis] - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: extent = np.array([0.12, 0.07, 0.07]) # m - return self._compile_scenario_1_problem(extent) + return self._compile_scenario_1_problem(extent, center_frequency) def get_default_source(self) -> sources.Source: """Return the default source for the scenario.""" @@ -252,11 +234,11 @@ def _create_scenario_1_mask(material, grid): elif material == "skin": _fill_mask(mask, start=interfaces[0], end=interfaces[1], dx=dx) - elif material == "cortical bone": + elif material == "cortical_bone": _fill_mask(mask, start=interfaces[1], end=interfaces[2], dx=dx) _fill_mask(mask, start=interfaces[3], end=interfaces[4], dx=dx) - elif material == "trabecular bone": + elif material == "trabecular_bone": _fill_mask(mask, start=interfaces[2], end=interfaces[3], dx=dx) elif material == "brain": diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py index a0ee1406..647f81d9 100644 --- a/src/neurotechdevkit/scenarios/_scenario_2.py +++ b/src/neurotechdevkit/scenarios/_scenario_2.py @@ -1,45 +1,21 @@ from __future__ import annotations import pathlib -from typing import Mapping, Protocol +from typing import Mapping import hdf5storage import numpy as np import numpy.typing as npt import stride -from mosaic.types import Struct from .. import rendering, sources -from . import materials -from ._base import Scenario2D, Scenario3D, Target +from ..materials import Material +from ._base import Scenario, Scenario2D, Scenario3D, Target from ._utils import add_material_fields_to_problem, make_grid -class _Scenario2MixinProtocol(Protocol): - """Provide type-hinting for Scenario 2 members used by mixins.""" - - @property - def scenario_id(self) -> str: - ... - - @property - def complexity(self) -> str: - ... - - @property - def materials(self) -> Mapping[str, Struct]: - ... - - @property - def layer_ids(self) -> Mapping[str, int]: - ... - - def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]: - ... - - -class _Scenario2Mixin: - """A mixin providing specific implementation detail for scenario 2. +class Scenario2(Scenario): + """Specific implementation detail for scenario 2. Scenario 2 is based on benchmark 8 of the following paper: @@ -51,30 +27,38 @@ class _Scenario2Mixin: https://asa.scitation.org/doi/pdf/10.1121/10.0013426 """ - @property - def _material_layers(self: _Scenario2MixinProtocol) -> list[tuple[str, Struct]]: - return [ - ("water", materials.water), - ("skull", materials.cortical_bone), - ("brain", materials.brain), - ] + material_layers = [ + "water", + "cortical_bone", + "brain", + ] + material_properties = { + "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"), + "cortical_bone": Material( + vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA" + ), + "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"), + } @property def _material_outline_upsample_factor(self) -> int: return 4 + def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]: + """Will be implemented by the subclass.""" + raise NotImplementedError() + def _compile_scenario_2_problem( - self: _Scenario2MixinProtocol, extent: npt.NDArray[np.float_] + self, extent: npt.NDArray[np.float_], center_frequency: float ) -> stride.Problem: # scenario constants speed_water = 1500 # m/s - c_freq = 500e3 # hz # desired resolution for complexity=fast ppw = 6 # compute resolution - dx = speed_water / c_freq / ppw # m + dx = speed_water / center_frequency / ppw # m grid = make_grid(extent=extent, dx=dx) problem = stride.Problem( @@ -82,14 +66,14 @@ def _compile_scenario_2_problem( ) problem = add_material_fields_to_problem( problem=problem, - materials=self.materials, + materials=self.get_materials(center_frequency), layer_ids=self.layer_ids, masks=self._get_material_masks(), ) return problem -class Scenario2_2D(_Scenario2Mixin, Scenario2D): +class Scenario2_2D(Scenario2, Scenario2D): """A 2D implementation of scenario 2. Scenario 2 is based on benchmark 8 of the following paper: @@ -148,14 +132,14 @@ def __init__(self, complexity="fast"): complexity=complexity, ) - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: extent = np.array([0.225, 0.170]) # m - return self._compile_scenario_2_problem(extent) + return self._compile_scenario_2_problem(extent, center_frequency) def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]: return { name: _create_scenario_2_mask(name, convert_2d=True) - for name in self.materials.keys() + for name in self.material_layers } def get_default_source(self) -> sources.Source: @@ -169,7 +153,7 @@ def get_default_source(self) -> sources.Source: ) -class Scenario2_3D(_Scenario2Mixin, Scenario3D): +class Scenario2_3D(Scenario2, Scenario3D): """A 3D implementation of scenario 2. Scenario 2 is based on benchmark 8 of the following paper: @@ -260,12 +244,12 @@ def viewer_config_3d(self) -> rendering.ViewerConfig3D: init_zoom=2.0, colormaps={ "water": "blue", - "skull": "magma", + "cortical_bone": "magma", "brain": "bop orange", }, opacities={ "water": 0.8, - "skull": 0.2, + "cortical_bone": 0.2, "brain": 0.2, }, ) @@ -279,14 +263,14 @@ def get_default_slice_position(self, axis: int) -> float: default_positions = np.array([0.1, 0.0, 0.0]) return default_positions[axis] - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: extent = np.array([0.225, 0.170, 0.190]) # m - return self._compile_scenario_2_problem(extent) + return self._compile_scenario_2_problem(extent, center_frequency) def _get_material_masks(self): return { name: _create_scenario_2_mask(name, convert_2d=False) - for name in self.materials.keys() + for name in self.material_layers } def get_default_source(self): @@ -309,7 +293,7 @@ def _create_scenario_2_mask(material, convert_2d=False) -> npt.NDArray[np.bool_] skull_mask = mat_data["skull_mask"].astype(np.bool_) brain_mask = mat_data["brain_mask"].astype(np.bool_) - if material == "skull": + if material == "cortical_bone": mask = skull_mask elif material == "brain": diff --git a/src/neurotechdevkit/scenarios/materials.py b/src/neurotechdevkit/scenarios/materials.py deleted file mode 100644 index 78222f25..00000000 --- a/src/neurotechdevkit/scenarios/materials.py +++ /dev/null @@ -1,48 +0,0 @@ -"""Materials for the neurotechdevkit scenarios.""" -from mosaic.types import Struct - -# TODO: encapsulate mosaic struct behind an NDK materials type - - -water = Struct() -water.vp = 1500.0 -water.rho = 1000.0 -water.alpha = 0.0 -water.render_color = "#2E86AB" - - -skin = Struct() -skin.vp = 1610.0 -skin.rho = 1090.0 -skin.alpha = 0.2 -skin.render_color = "#FA8B53" - - -cortical_bone = Struct() -cortical_bone.vp = 2800.0 -cortical_bone.rho = 1850.0 -cortical_bone.alpha = 4.0 -cortical_bone.render_color = "#FAF0CA" - - -trabecular_bone = Struct() -trabecular_bone.vp = 2300.0 -trabecular_bone.rho = 1700.0 -trabecular_bone.alpha = 8.0 -trabecular_bone.render_color = "#EBD378" - - -brain = Struct() -brain.vp = 1560.0 -brain.rho = 1040.0 -brain.alpha = 0.3 -brain.render_color = "#DB504A" - - -# these numbers are completely made up -# TODO: research reasonable values -tumor = Struct() -tumor.vp = 1650.0 -tumor.rho = 1150.0 -tumor.alpha = 0.8 -tumor.render_color = "#94332F" diff --git a/tests/neurotechdevkit/scenarios/test_base.py b/tests/neurotechdevkit/scenarios/test_base.py index c4454024..fc2af5a9 100644 --- a/tests/neurotechdevkit/scenarios/test_base.py +++ b/tests/neurotechdevkit/scenarios/test_base.py @@ -5,10 +5,8 @@ import pytest import stride from frozenlist import FrozenList -from mosaic.types import Struct from neurotechdevkit.results import PulsedResult, SteadyStateResult -from neurotechdevkit.scenarios import materials from neurotechdevkit.scenarios._base import Scenario, Target from neurotechdevkit.scenarios._utils import make_grid, wavelet_helper from neurotechdevkit.sources import FocusedSource3D, PlanarSource3D, Source @@ -32,6 +30,7 @@ class ScenarioBaseTester(Scenario): description="bar", ), } + material_layers = ["brain", "skin"] default_source = PlanarSource3D( position=np.array([0.05, 0.1, 0.05]), @@ -42,20 +41,17 @@ class ScenarioBaseTester(Scenario): def __init__(self, complexity="fast"): self._target_id = "target_1" + self._problem = self._compile_problem(center_frequency=5e5) super().__init__( origin=np.array([-0.1, -0.1, 0.0]), complexity=complexity, ) - @property - def _material_layers(self) -> list[tuple[str, Struct]]: - return [("foo", materials.brain), ("bar", materials.skin)] - @property def _material_outline_upsample_factor(self) -> int: return 3 - def _compile_problem(self) -> stride.Problem: + def _compile_problem(self, center_frequency: float) -> stride.Problem: extent = np.array([2.0, 3.0, 4.0]) dx = 0.1 grid = make_grid(extent=extent, dx=dx) @@ -519,7 +515,7 @@ def test_get_layer_mask_with_wrong_layer_name(tester_with_layers): def test_get_layer_mask_with_first_layer(tester_with_layers): """Verify that get_layer_mask returns the expected mask for the first layer.""" - mask = tester_with_layers.get_layer_mask("foo") + mask = tester_with_layers.get_layer_mask("brain") expected = np.zeros_like(mask, dtype=bool) expected[:5] = True np.testing.assert_allclose(mask, expected) @@ -527,7 +523,7 @@ def test_get_layer_mask_with_first_layer(tester_with_layers): def test_get_layer_mask_with_last_layer(tester_with_layers): """Verify that get_layer_mask returns the expected mask for the last layer.""" - mask = tester_with_layers.get_layer_mask("bar") + mask = tester_with_layers.get_layer_mask("skin") expected = np.zeros_like(mask, dtype=bool) expected[5:] = True np.testing.assert_allclose(mask, expected) diff --git a/tests/neurotechdevkit/scenarios/test_materials.py b/tests/neurotechdevkit/scenarios/test_materials.py new file mode 100644 index 00000000..e4728a6d --- /dev/null +++ b/tests/neurotechdevkit/scenarios/test_materials.py @@ -0,0 +1,108 @@ +import numpy as np +import pytest +from mosaic.types import Struct + +from neurotechdevkit.materials import Material +from neurotechdevkit.scenarios import Scenario + + +def compare_structs(struct1: Struct, struct2: Struct): + """Compare two Structs.""" + assert struct1.vp == struct2.vp + assert struct1.rho == struct2.rho + assert struct1.alpha == struct2.alpha + assert struct1.render_color == struct2.render_color + + +class BaseScenario(Scenario): + """A scenario for testing the materials module.""" + + _SCENARIO_ID = "scenario-tester" + _TARGET_OPTIONS = {} + + def __init__(self, complexity="fast"): + self._target_id = "target_1" + super().__init__( + origin=np.array([-0.1, -0.1, 0.0]), + complexity=complexity, + ) + + def _compile_problem(self, center_frequency: float): + pass + + def _material_outline_upsample_factor(self): + pass + + def get_default_source(self): + pass + + def get_target_mask(self): + pass + + +def test_custom_material_property(): + """Test that a custom material property is used.""" + + class ScenarioWithCustomMaterialProperties(BaseScenario): + material_layers = ["brain"] + material_properties = { + "brain": Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB") + } + + scenario = ScenarioWithCustomMaterialProperties() + assert scenario.material_colors == {"brain": "#2E86AB"} + + materials = scenario.get_materials(500e3) + assert list(materials.keys()) == ["brain"] + compare_structs( + materials["brain"], + Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(), + ) + + +def test_new_material(): + """Test that a new material is used.""" + + class ScenarioWithCustomMaterial(BaseScenario): + material_layers = ["brain", "eye"] + material_properties = { + "eye": Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB") + } + + scenario = ScenarioWithCustomMaterial() + assert scenario.material_colors == {"brain": "#DB504A", "eye": "#2E86AB"} + + materials = scenario.get_materials(500e3) + + assert list(materials.keys()) == ["brain", "eye"] + compare_structs( + materials["eye"], + Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(), + ) + + +def test_material_absorption_is_calculated(): + """Test that the material absorption is calculated for a frequency !=500e3.""" + + class ScenarioWithBrainMaterial(BaseScenario): + material_layers = ["brain"] + + scenario = ScenarioWithBrainMaterial() + + materials = scenario.get_materials(600e3) + + assert list(materials.keys()) == ["brain"] + assert materials["brain"].alpha == 0.3041793231753331 + + +def test_unknown_material_without_properties(): + """Test that an unknown material without properties raises an error.""" + + class ScenarioWithCustomMaterial(BaseScenario): + material_layers = ["unknown_material"] + + scenario = ScenarioWithCustomMaterial() + with pytest.raises(ValueError): + scenario.material_colors + with pytest.raises(ValueError): + scenario.get_materials(500e3) diff --git a/tests/neurotechdevkit/scenarios/test_metrics.py b/tests/neurotechdevkit/scenarios/test_metrics.py index d6f56fef..baaf91fb 100644 --- a/tests/neurotechdevkit/scenarios/test_metrics.py +++ b/tests/neurotechdevkit/scenarios/test_metrics.py @@ -54,7 +54,7 @@ def fake_result(fake_steady_state_matrix): """Returns a fake SteadyStateResult2D that can be used for testing.""" return SteadyStateResult2D( scenario=SimpleNamespace( - ordered_layers=["water", "brain"], + material_layers=["water", "brain"], get_layer_mask=fake_layer_masks, get_target_mask=lambda: fake_layer_masks("target"), get_field_data=fake_fields, diff --git a/tests/neurotechdevkit/scenarios/test_time.py b/tests/neurotechdevkit/scenarios/test_time.py index 30c4133a..d9f2d4cf 100644 --- a/tests/neurotechdevkit/scenarios/test_time.py +++ b/tests/neurotechdevkit/scenarios/test_time.py @@ -52,7 +52,7 @@ def test_select_simulation_time_for_steady_state_with_defined_time_to_ss(grid_fa When time_to_steady_state is provided, the function should use it to calculate the simulation time. """ - test_materials = {"test": materials.brain} + test_materials = {"brain": materials.get_material("brain")} test_freq = 2.0e4 n_cycles = 5 test_time_ss = 1e-4 @@ -76,7 +76,10 @@ def test_select_simulation_time_for_steady_state_with_default_time_to_ss(grid_fa When time_to_steady_state is not provided, the function should estimate a value for it based on the size of the scenario and the lowest speed of sound. """ - test_materials = {"fast": materials.cortical_bone, "slow": materials.brain} + test_materials = { + "fast": materials.get_material("cortical_bone"), + "slow": materials.get_material("brain"), + } test_freq = 2.0e4 n_cycles = 5 sim_time = select_simulation_time_for_steady_state( @@ -89,14 +92,17 @@ def test_select_simulation_time_for_steady_state_with_default_time_to_ss(grid_fa ) period = 1.0 / test_freq length = np.sqrt(0.3**2 + 0.4**2 + 0.5**2) - expected_time_ss = 2 * length / materials.brain.vp + expected_time_ss = 2 * length / materials.get_material("brain").vp expected_time = expected_time_ss + n_cycles * period np.testing.assert_allclose(sim_time, expected_time) def test_select_simulation_time_for_pulsed(grid_fake): """Verify that the returned simulation time is correct for a pulsed simulation.""" - test_materials = {"fast": materials.cortical_bone, "slow": materials.brain} + test_materials = { + "fast": materials.get_material("cortical_bone"), + "slow": materials.get_material("brain"), + } source_delay = 31.0 sim_time = select_simulation_time_for_pulsed( grid=grid_fake, @@ -104,7 +110,7 @@ def test_select_simulation_time_for_pulsed(grid_fake): delay=source_delay, ) length = np.sqrt(0.3**2 + 0.4**2 + 0.5**2) - expected_time = source_delay + length / materials.brain.vp + expected_time = source_delay + length / materials.get_material("brain").vp np.testing.assert_allclose(sim_time, expected_time) diff --git a/whitelist.txt b/whitelist.txt index 899953c8..84f3b8f6 100644 --- a/whitelist.txt +++ b/whitelist.txt @@ -1,3 +1,4 @@ +α al antiparallel args