diff --git a/pytissueoptics/examples/benchmarks/cube60.py b/pytissueoptics/examples/benchmarks/cube60.py index 9ad0ef7e..1ae095c7 100644 --- a/pytissueoptics/examples/benchmarks/cube60.py +++ b/pytissueoptics/examples/benchmarks/cube60.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "MCX Homogeneous cube" diff --git a/pytissueoptics/examples/benchmarks/cubesph60b.py b/pytissueoptics/examples/benchmarks/cubesph60b.py index 7f30f8c0..2f490ef8 100644 --- a/pytissueoptics/examples/benchmarks/cubesph60b.py +++ b/pytissueoptics/examples/benchmarks/cubesph60b.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "MCX Sphere" diff --git a/pytissueoptics/examples/benchmarks/skinvessel.py b/pytissueoptics/examples/benchmarks/skinvessel.py index efc57603..2829790c 100644 --- a/pytissueoptics/examples/benchmarks/skinvessel.py +++ b/pytissueoptics/examples/benchmarks/skinvessel.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "Skin vessel" diff --git a/pytissueoptics/examples/benchmarks/sphshells.py b/pytissueoptics/examples/benchmarks/sphshells.py index 9e7ae6b9..d912f817 100644 --- a/pytissueoptics/examples/benchmarks/sphshells.py +++ b/pytissueoptics/examples/benchmarks/sphshells.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "MCX Spherical shells" diff --git a/pytissueoptics/examples/rayscattering/ex01.py b/pytissueoptics/examples/rayscattering/ex01.py index 2b142023..7f85c290 100644 --- a/pytissueoptics/examples/rayscattering/ex01.py +++ b/pytissueoptics/examples/rayscattering/ex01.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "Divergent source propagation through a multi-layered tissue" @@ -24,6 +25,7 @@ def exampleCode(): viewer.reportStats() viewer.show2D(View2DProjectionX()) + viewer.show2D(View2DProjectionX(energyType=EnergyType.FLUENCE_RATE)) viewer.show2D(View2DProjectionX(solidLabel="middleLayer")) viewer.show2D(View2DSurfaceZ(solidLabel="middleLayer", surfaceLabel="interface1", surfaceEnergyLeaving=False)) viewer.show1D(Direction.Z_POS) diff --git a/pytissueoptics/examples/rayscattering/ex02.py b/pytissueoptics/examples/rayscattering/ex02.py index d849a8a6..a5299ddf 100644 --- a/pytissueoptics/examples/rayscattering/ex02.py +++ b/pytissueoptics/examples/rayscattering/ex02.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "Propagation in an Infinite Medium" diff --git a/pytissueoptics/examples/rayscattering/ex03.py b/pytissueoptics/examples/rayscattering/ex03.py index 18a79aee..0ccceac8 100644 --- a/pytissueoptics/examples/rayscattering/ex03.py +++ b/pytissueoptics/examples/rayscattering/ex03.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = ( diff --git a/pytissueoptics/examples/rayscattering/ex04.py b/pytissueoptics/examples/rayscattering/ex04.py index 532e5efe..3720e09c 100644 --- a/pytissueoptics/examples/rayscattering/ex04.py +++ b/pytissueoptics/examples/rayscattering/ex04.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "Custom layer stack" diff --git a/pytissueoptics/examples/rayscattering/ex05.py b/pytissueoptics/examples/rayscattering/ex05.py index 450a6110..b798cb0c 100644 --- a/pytissueoptics/examples/rayscattering/ex05.py +++ b/pytissueoptics/examples/rayscattering/ex05.py @@ -1,4 +1,5 @@ import env # noqa: F401 + from pytissueoptics import * # noqa: F403 TITLE = "Sphere inside a cube" diff --git a/pytissueoptics/rayscattering/__init__.py b/pytissueoptics/rayscattering/__init__.py index 8ab8c01e..4abdcb29 100644 --- a/pytissueoptics/rayscattering/__init__.py +++ b/pytissueoptics/rayscattering/__init__.py @@ -13,7 +13,7 @@ View2DSurfaceY, View2DSurfaceZ, ) -from .energyLogging import EnergyLogger +from .energyLogging import EnergyLogger, EnergyType from .materials import ScatteringMaterial from .opencl import CONFIG, disableOpenCL, hardwareAccelerationIsAvailable from .photon import Photon @@ -29,6 +29,7 @@ "DirectionalSource", "DivergentSource", "EnergyLogger", + "EnergyType", "ScatteringScene", "Viewer", "PointCloudStyle", diff --git a/pytissueoptics/rayscattering/display/profiles/profile1D.py b/pytissueoptics/rayscattering/display/profiles/profile1D.py index 67a4ca9d..19a4fac4 100644 --- a/pytissueoptics/rayscattering/display/profiles/profile1D.py +++ b/pytissueoptics/rayscattering/display/profiles/profile1D.py @@ -4,20 +4,29 @@ from matplotlib import pyplot as plt from pytissueoptics.rayscattering.display.utils import Direction +from pytissueoptics.rayscattering.energyLogging import EnergyType class Profile1D: """ Since 1D profiles are easily generated from existing 2D views or 3D data, this class is only used as a small dataclass. Only used internally Profile1DFactory when Viewer.show1D() is called. The user should only use the - endpoint Viewer.show1D() which doesn't require to create a Profile1D object. + endpoint Viewer.show1D() which doesn't require creating a Profile1D object. """ - def __init__(self, data: np.ndarray, horizontalDirection: Direction, limits: Tuple[float, float], name: str = None): + def __init__( + self, + data: np.ndarray, + horizontalDirection: Direction, + limits: Tuple[float, float], + name: str = None, + energyType=EnergyType.DEPOSITION, + ): self.data = data self.limits = limits self.horizontalDirection = horizontalDirection self.name = name + self.energyType = energyType def show(self, logScale: bool = True): limits = sorted(self.limits) @@ -33,5 +42,5 @@ def show(self, logScale: bool = True): plt.title(self.name) plt.xlim(*limits) plt.xlabel("xyz"[self.horizontalDirection.axis]) - plt.ylabel("Energy") + plt.ylabel("Deposited energy" if self.energyType == EnergyType.DEPOSITION else "Fluence rate") plt.show() diff --git a/pytissueoptics/rayscattering/display/profiles/profileFactory.py b/pytissueoptics/rayscattering/display/profiles/profileFactory.py index 765276a3..6a8d1ad6 100644 --- a/pytissueoptics/rayscattering/display/profiles/profileFactory.py +++ b/pytissueoptics/rayscattering/display/profiles/profileFactory.py @@ -5,7 +5,7 @@ from pytissueoptics.rayscattering import utils from pytissueoptics.rayscattering.display.profiles import Profile1D from pytissueoptics.rayscattering.display.utils import Direction -from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloudFactory +from pytissueoptics.rayscattering.energyLogging import EnergyLogger, EnergyType, PointCloudFactory from pytissueoptics.rayscattering.scatteringScene import ScatteringScene @@ -28,6 +28,7 @@ def create( surfaceEnergyLeaving: bool = True, limits: Tuple[float, float] = None, binSize: float = None, + energyType=EnergyType.DEPOSITION, ) -> Profile1D: solidLabel, surfaceLabel = self._correctCapitalization(solidLabel, surfaceLabel) if binSize is None: @@ -39,15 +40,15 @@ def create( if self._logger.has3D: histogram = self._extractHistogramFrom3D( - horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, bins + horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, bins, energyType ) else: histogram = self._extractHistogramFromViews( - horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, bins + horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, bins, energyType ) - name = self._createName(horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving) - return Profile1D(histogram, horizontalDirection, limits, name) + name = self._createName(horizontalDirection, solidLabel, surfaceLabel, surfaceEnergyLeaving, energyType) + return Profile1D(histogram, horizontalDirection, limits, name, energyType) def _getDefaultLimits(self, horizontalDirection: Direction, solidLabel: str = None): if solidLabel: @@ -70,8 +71,9 @@ def _extractHistogramFrom3D( surfaceEnergyLeaving: bool, limits: Tuple[float, float], bins: int, + energyType: EnergyType, ): - pointCloud = self._pointCloudFactory.getPointCloud(solidLabel, surfaceLabel) + pointCloud = self._pointCloudFactory.getPointCloud(solidLabel, surfaceLabel, energyType=energyType) if surfaceLabel: if surfaceEnergyLeaving: @@ -95,7 +97,10 @@ def _extractHistogramFromViews( surfaceEnergyLeaving: bool, limits: Tuple[float, float], bins: int, + energyType: EnergyType, ): + energyMismatch = False + for view in self._logger.views: if view.axis == horizontalDirection.axis: continue @@ -118,10 +123,16 @@ def _extractHistogramFromViews( continue if viewBins != bins: continue + if not surfaceLabel and energyType != view.energyType: + energyMismatch = True + continue return self._extractHistogramFromView(view, horizontalDirection) - raise RuntimeError("Cannot create 1D profile. The 3D data was discarded and no matching 2D view was found.") + error_message = "Cannot create 1D profile. The 3D data was discarded and no matching 2D view was found." + if energyMismatch: + error_message += " Note that a view candidate was found to only differ in energy type." + raise RuntimeError(error_message) def _extractHistogramFromView(self, view, horizontalDirection: Direction): if view.axisU == horizontalDirection.axis: @@ -160,9 +171,14 @@ def _correctCapitalization(self, solidLabel, surfaceLabel): return solidLabel, surfaceLabel def _createName( - self, horizontalDirection: Direction, solidLabel: str, surfaceLabel: str, surfaceEnergyLeaving: bool + self, + horizontalDirection: Direction, + solidLabel: str, + surfaceLabel: str, + surfaceEnergyLeaving: bool, + energyType: EnergyType, ) -> str: - name = "Energy profile along " + "xyz"[horizontalDirection.axis] + name = f"{energyType.name} profile along " + "xyz"[horizontalDirection.axis] if solidLabel: name += " of " + solidLabel if surfaceLabel: diff --git a/pytissueoptics/rayscattering/display/viewer.py b/pytissueoptics/rayscattering/display/viewer.py index 614ca2e4..adf529cb 100644 --- a/pytissueoptics/rayscattering/display/viewer.py +++ b/pytissueoptics/rayscattering/display/viewer.py @@ -7,7 +7,11 @@ from pytissueoptics.rayscattering.display.profiles import ProfileFactory from pytissueoptics.rayscattering.display.utils import Direction from pytissueoptics.rayscattering.display.views import View2D, ViewGroup -from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloudFactory +from pytissueoptics.rayscattering.energyLogging import ( + EnergyLogger, + EnergyType, + PointCloudFactory, +) from pytissueoptics.rayscattering.energyLogging.pointCloud import PointCloud from pytissueoptics.rayscattering.scatteringScene import ScatteringScene from pytissueoptics.rayscattering.source import Source @@ -40,6 +44,7 @@ class PointCloudStyle: showSolidPoints (bool): Show the point clouds of the solids. showSurfacePointsLeaving (bool): Show energy that left the surface (direction with surface normal). showSurfacePointsEntering (bool): Show energy that entered the surface (direction opposite to surface normal). + energyType (EnergyType): Type of energy to show for volumetric datapoints (deposition or fluence). Other attributes: showPointsAsSpheres (bool): Show the points as spheres or as dots. Dots require less memory. @@ -60,6 +65,7 @@ def __init__( showSolidPoints: bool = True, showSurfacePointsLeaving: bool = True, showSurfacePointsEntering: bool = False, + energyType=EnergyType.DEPOSITION, showPointsAsSpheres: bool = False, pointSize: float = 0.15, scaleWithValue: bool = True, @@ -75,6 +81,7 @@ def __init__( self.showSolidPoints = showSolidPoints self.showSurfacePointsLeaving = showSurfacePointsLeaving self.showSurfacePointsEntering = showSurfacePointsEntering + self.energyType = energyType self.showPointsAsSpheres = showPointsAsSpheres self.pointSize = pointSize @@ -142,6 +149,7 @@ def show3DVolumeSlicer( logScale: bool = True, interpolate: bool = False, limits: Tuple[tuple, tuple, tuple] = None, + energyType: EnergyType = EnergyType.DEPOSITION, ): if not self._logger.has3D: utils.warn("ERROR: Cannot show 3D volume slicer without 3D data.") @@ -161,7 +169,7 @@ def show3DVolumeSlicer( f"Consider using a larger binSize or tighter limits." ) - points = self._pointCloudFactory.getPointCloudOfSolids().solidPoints + points = self._pointCloudFactory.getPointCloudOfSolids(energyType=energyType).solidPoints try: hist, _ = np.histogramdd(points[:, 1:], bins=bins, weights=points[:, 0], range=limits) except MemoryError: @@ -195,8 +203,11 @@ def show1D( surfaceEnergyLeaving: bool = True, limits: Tuple[float, float] = None, binSize: float = None, + energyType: EnergyType = EnergyType.DEPOSITION, ): - profile = self._profileFactory.create(along, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, binSize) + profile = self._profileFactory.create( + along, solidLabel, surfaceLabel, surfaceEnergyLeaving, limits, binSize, energyType + ) profile.show(logScale=logScale) def reportStats(self, solidLabel: str = None, saveToFile: str = None, verbose=True): @@ -204,7 +215,9 @@ def reportStats(self, solidLabel: str = None, saveToFile: str = None, verbose=Tr stats.report(solidLabel=solidLabel, saveToFile=saveToFile, verbose=verbose) def _addPointCloud(self, style: PointCloudStyle): - pointCloud = self._pointCloudFactory.getPointCloud(solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel) + pointCloud = self._pointCloudFactory.getPointCloud( + solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel, energyType=style.energyType + ) self._drawPointCloudOfSolids(pointCloud, style) self._drawPointCloudOfSurfaces(pointCloud, style) diff --git a/pytissueoptics/rayscattering/display/views/defaultViews.py b/pytissueoptics/rayscattering/display/views/defaultViews.py index 1f316811..a18092eb 100644 --- a/pytissueoptics/rayscattering/display/views/defaultViews.py +++ b/pytissueoptics/rayscattering/display/views/defaultViews.py @@ -9,6 +9,7 @@ Direction, ) from pytissueoptics.rayscattering.display.views.view2D import View2D, ViewGroup +from pytissueoptics.rayscattering.energyLogging.energyType import EnergyType class View2DProjection(View2D): @@ -19,9 +20,15 @@ def __init__( solidLabel: str = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): super().__init__( - projectionDirection, horizontalDirection, solidLabel=solidLabel, limits=limits, binSize=binSize + projectionDirection, + horizontalDirection, + solidLabel=solidLabel, + limits=limits, + binSize=binSize, + energyType=energyType, ) def _filter(self, dataPoints: np.ndarray) -> np.ndarray: @@ -34,8 +41,11 @@ def __init__( solidLabel: str = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): - super().__init__(*DEFAULT_X_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize) + super().__init__( + *DEFAULT_X_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + ) class View2DProjectionY(View2DProjection): @@ -44,8 +54,11 @@ def __init__( solidLabel: str = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): - super().__init__(*DEFAULT_Y_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize) + super().__init__( + *DEFAULT_Y_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + ) class View2DProjectionZ(View2DProjection): @@ -54,8 +67,11 @@ def __init__( solidLabel: str = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): - super().__init__(*DEFAULT_Z_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize) + super().__init__( + *DEFAULT_Z_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + ) class View2DSurface(View2D): @@ -159,6 +175,7 @@ def __init__( thickness: float = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): super().__init__( projectionDirection, @@ -168,6 +185,7 @@ def __init__( thickness=thickness, limits=limits, binSize=binSize, + energyType=energyType, ) def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float, float, float]): @@ -191,6 +209,7 @@ def __init__( thickness: float = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -199,6 +218,7 @@ def __init__( thickness=thickness, limits=limits, binSize=binSize, + energyType=energyType, ) @@ -210,6 +230,7 @@ def __init__( thickness: float = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -218,6 +239,7 @@ def __init__( thickness=thickness, limits=limits, binSize=binSize, + energyType=energyType, ) @@ -229,6 +251,7 @@ def __init__( thickness: float = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, @@ -237,4 +260,5 @@ def __init__( thickness=thickness, limits=limits, binSize=binSize, + energyType=energyType, ) diff --git a/pytissueoptics/rayscattering/display/views/view2D.py b/pytissueoptics/rayscattering/display/views/view2D.py index 021edb13..e61e6de2 100644 --- a/pytissueoptics/rayscattering/display/views/view2D.py +++ b/pytissueoptics/rayscattering/display/views/view2D.py @@ -13,6 +13,7 @@ DEFAULT_Z_VIEW_DIRECTIONS, Direction, ) +from pytissueoptics.rayscattering.energyLogging.energyType import EnergyType class ViewGroup(Flag): @@ -48,6 +49,7 @@ def __init__( thickness: float = None, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + energyType=EnergyType.DEPOSITION, ): """ The 2D view plane is obtained by looking towards the 'projectionDirection'. The 'horizontalDirection' @@ -75,6 +77,7 @@ def __init__( self._surfaceEnergyLeaving = surfaceEnergyLeaving self._position = position self._thickness = thickness + self._energyType = energyType limits = [sorted(limit) for limit in limits] if limits else [None, None] self._limitsU, self._limitsV = limits @@ -114,6 +117,10 @@ def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float "Cannot allocate memory for 2D view. Consider increasing `defaultBinSize` of EnergyLogger." ) + @property + def energyType(self) -> EnergyType: + return self._energyType + def extractData(self, dataPoints: np.ndarray): """ Used internally by Logger2D to store 3D datapoints into this 2D view. @@ -236,6 +243,8 @@ def isEqualTo(self, other: "View2D") -> bool: otherLimits = sorted(other._limitsU), sorted(other._limitsV) if limits != otherLimits: return False + if self._energyType != other._energyType: + return False return True def isContainedBy(self, other: "View2D") -> bool: @@ -346,7 +355,8 @@ def name(self) -> str: if self._surfaceLabel: objectLabel += f" surface {self._surfaceLabel}" objectLabel += " (leaving)" if self._surfaceEnergyLeaving else " (entering)" - + elif self._energyType == EnergyType.FLUENCE_RATE: + objectLabel += " (fluence)" return f"{self.__class__.__name__} of {objectLabel}" @property diff --git a/pytissueoptics/rayscattering/display/views/viewFactory.py b/pytissueoptics/rayscattering/display/views/viewFactory.py index 01113c97..edd00753 100644 --- a/pytissueoptics/rayscattering/display/views/viewFactory.py +++ b/pytissueoptics/rayscattering/display/views/viewFactory.py @@ -12,14 +12,20 @@ View2DSurfaceZ, ) from pytissueoptics.rayscattering.display.views.view2D import View2D, ViewGroup +from pytissueoptics.rayscattering.energyLogging.energyType import EnergyType from pytissueoptics.rayscattering.scatteringScene import ScatteringScene class ViewFactory: def __init__( - self, scene: ScatteringScene, defaultBinSize: Union[float, Tuple[float, float, float]], infiniteLimits: tuple + self, + scene: ScatteringScene, + defaultBinSize: Union[float, Tuple[float, float, float]], + infiniteLimits: tuple, + energyType=EnergyType.DEPOSITION, ): self._scene = scene + self._energyType = energyType self._defaultBinSize3D = defaultBinSize if isinstance(self._defaultBinSize3D, float): @@ -137,10 +143,9 @@ def _viewHasValidSurfaceLabel(self, view) -> bool: return self._viewHasValidSurfaceLabel(view) return False - @staticmethod - def _getDefaultViewsXYZ(solidLabel: str = None) -> List[View2D]: + def _getDefaultViewsXYZ(self, solidLabel: str = None) -> List[View2D]: return [ - View2DProjectionX(solidLabel=solidLabel), - View2DProjectionY(solidLabel=solidLabel), - View2DProjectionZ(solidLabel=solidLabel), + View2DProjectionX(solidLabel=solidLabel, energyType=self._energyType), + View2DProjectionY(solidLabel=solidLabel, energyType=self._energyType), + View2DProjectionZ(solidLabel=solidLabel, energyType=self._energyType), ] diff --git a/pytissueoptics/rayscattering/energyLogging/__init__.py b/pytissueoptics/rayscattering/energyLogging/__init__.py index a36c755d..c9048460 100644 --- a/pytissueoptics/rayscattering/energyLogging/__init__.py +++ b/pytissueoptics/rayscattering/energyLogging/__init__.py @@ -1,9 +1,11 @@ from .energyLogger import EnergyLogger +from .energyType import EnergyType from .pointCloud import PointCloud from .pointCloudFactory import PointCloudFactory __all__ = [ "EnergyLogger", + "EnergyType", "PointCloud", "PointCloudFactory", ] diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index ed63576b..d4be2558 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -1,6 +1,6 @@ import os import pickle -from typing import List, Union +from typing import List, Union, Optional import numpy as np @@ -9,7 +9,10 @@ from pytissueoptics.rayscattering.display.views.viewFactory import ViewFactory from pytissueoptics.rayscattering.scatteringScene import ScatteringScene from pytissueoptics.scene.geometry import Vector -from pytissueoptics.scene.logger.logger import InteractionKey, Logger +from pytissueoptics.scene.logger.listArrayContainer import ListArrayContainer +from pytissueoptics.scene.logger.logger import DataType, InteractionKey, Logger + +from .energyType import EnergyType class EnergyLogger(Logger): @@ -19,6 +22,7 @@ def __init__( filepath: str = None, keep3D: bool = True, views: Union[ViewGroup, List[View2D]] = ViewGroup.ALL, + defaultViewEnergyType: EnergyType = EnergyType.DEPOSITION, defaultBinSize: Union[float, tuple] = 0.01, infiniteLimits=((-5, 5), (-5, 5), (-5, 5)), ): @@ -47,7 +51,7 @@ def __init__( self._keep3D = keep3D self._defaultBinSize = defaultBinSize self._infiniteLimits = infiniteLimits - self._viewFactory = ViewFactory(scene, defaultBinSize, infiniteLimits) + self._viewFactory = ViewFactory(scene, defaultBinSize, infiniteLimits, energyType=defaultViewEnergyType) self._sceneHash = hash(scene) self._defaultViews = views @@ -193,6 +197,7 @@ def _getViewIndex(self, view: View2D) -> int: for i, v in enumerate(self._views): if v.isEqualTo(view): return i + raise ValueError(f"View {view.name} not found in the list of views.") def _viewExists(self, view: View2D) -> bool: return any([view.isEqualTo(v) for v in self._views]) @@ -229,7 +234,7 @@ def logDataPoint(self, value: float, position: Vector, key: InteractionKey): def _compileViews(self, views: List[View2D]): for key, data in self._data.items(): - datapointsContainer = data.dataPoints + datapointsContainer: Optional[ListArrayContainer] = data.dataPoints if datapointsContainer is None or len(datapointsContainer) == 0: continue for view in views: @@ -239,7 +244,12 @@ def _compileViews(self, views: List[View2D]): continue if view.surfaceLabel is None and key.surfaceLabel is not None: continue - view.extractData(datapointsContainer.getData()) + + data = datapointsContainer.getData() + if view.energyType == EnergyType.FLUENCE_RATE: + data = self._fluenceTransform(key, data) + + view.extractData(data) for view in views: self._outdatedViews.discard(view) @@ -272,3 +282,28 @@ def logSegment(self, start: Vector, end: Vector, key: InteractionKey = None): def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None): raise NotImplementedError("Can only log data points to an EnergyLogger.") + + def getDataPoints(self, key: InteractionKey, energyType=EnergyType.DEPOSITION) -> np.ndarray: + """All 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) where + the second axis is (value, x, y, z). The value can be the energy deposited, the fluence rate, or the + energy that crossed the surface. + + :param key: Filtering the data by solidLabel and surfaceLabel. + :param energyType: The type of volumetric energy to return when no surfaceLabel is given. + + :return: The data points (value, x, y, z) for the given solidLabel and surfaceLabel. If a surfaceLabel is given, + the value corresponds to the energy that crossed the surface (positive when in the direction of the normal). If + only a solidLabel is given, the value corresponds to the volumetric EnergyType at that point. + """ + if energyType == EnergyType.FLUENCE_RATE: + return self._getData(DataType.DATA_POINT, key, transform=self._fluenceTransform) + + return self._getData(DataType.DATA_POINT, key) + + def _fluenceTransform(self, key: InteractionKey, data: Optional[np.ndarray]) -> Optional[np.ndarray]: + # Converts volumetric data to fluence rate when needed. + if not key.volumetric or data is None: + return data + + data[:, 0] = data[:, 0] / self._scene.getMaterial(key.solidLabel).mu_a + return data diff --git a/pytissueoptics/rayscattering/energyLogging/energyType.py b/pytissueoptics/rayscattering/energyLogging/energyType.py new file mode 100644 index 00000000..9cb20154 --- /dev/null +++ b/pytissueoptics/rayscattering/energyLogging/energyType.py @@ -0,0 +1,10 @@ +from enum import Enum, auto + + +class EnergyType(Enum): + """ + Type of volumetric energy: either as the deposited energy in the solid (absorption) or as the fluence rate. + """ + + DEPOSITION = auto() + FLUENCE_RATE = auto() diff --git a/pytissueoptics/rayscattering/energyLogging/pointCloudFactory.py b/pytissueoptics/rayscattering/energyLogging/pointCloudFactory.py index 05486585..27cc48d9 100644 --- a/pytissueoptics/rayscattering/energyLogging/pointCloudFactory.py +++ b/pytissueoptics/rayscattering/energyLogging/pointCloudFactory.py @@ -1,25 +1,32 @@ import numpy as np -from pytissueoptics.rayscattering.energyLogging import PointCloud -from pytissueoptics.scene.logger import InteractionKey, Logger +from pytissueoptics.scene.logger import InteractionKey + +from .energyLogger import EnergyLogger +from .energyType import EnergyType +from .pointCloud import PointCloud class PointCloudFactory: - def __init__(self, logger: Logger): + def __init__(self, logger: EnergyLogger): self._logger = logger - def getPointCloud(self, solidLabel: str = None, surfaceLabel: str = None) -> PointCloud: + def getPointCloud( + self, solidLabel: str = None, surfaceLabel: str = None, energyType=EnergyType.DEPOSITION + ) -> PointCloud: if not solidLabel and not surfaceLabel: - return PointCloud(self.getPointCloudOfSolids().solidPoints, self.getPointCloudOfSurfaces().surfacePoints) - points = self._logger.getDataPoints(InteractionKey(solidLabel, surfaceLabel)) + return PointCloud( + self.getPointCloudOfSolids(energyType).solidPoints, self.getPointCloudOfSurfaces().surfacePoints + ) + points = self._logger.getDataPoints(InteractionKey(solidLabel, surfaceLabel), energyType=energyType) if surfaceLabel: return PointCloud(None, points) return PointCloud(points, None) - def getPointCloudOfSolids(self) -> PointCloud: + def getPointCloudOfSolids(self, energyType=EnergyType.DEPOSITION) -> PointCloud: points = [] for solidLabel in self._logger.getStoredSolidLabels(): - solidPoints = self.getPointCloud(solidLabel).solidPoints + solidPoints = self.getPointCloud(solidLabel, energyType=energyType).solidPoints if solidPoints is not None: points.append(solidPoints) if len(points) == 0: diff --git a/pytissueoptics/rayscattering/samples/phantomTissue.py b/pytissueoptics/rayscattering/samples/phantomTissue.py index fdf0a48d..3da47625 100644 --- a/pytissueoptics/rayscattering/samples/phantomTissue.py +++ b/pytissueoptics/rayscattering/samples/phantomTissue.py @@ -15,15 +15,15 @@ def __init__(self, worldMaterial=ScatteringMaterial()): def _create(self): n = [1.4, 1.7, 1.4] mu_s = [2, 3, 2] - mu_a = 1 + mu_a = [1, 1, 2] g = 0.8 w = 3 t = [0.75, 0.5, 0.75] - frontLayer = Cuboid(w, w, t[0], material=ScatteringMaterial(mu_s[0], mu_a, g, n[0]), label="frontLayer") - middleLayer = Cuboid(w, w, t[1], material=ScatteringMaterial(mu_s[1], mu_a, g, n[1]), label="middleLayer") - backLayer = Cuboid(w, w, t[2], material=ScatteringMaterial(mu_s[2], mu_a, g, n[2]), label="backLayer") + frontLayer = Cuboid(w, w, t[0], material=ScatteringMaterial(mu_s[0], mu_a[0], g, n[0]), label="frontLayer") + middleLayer = Cuboid(w, w, t[1], material=ScatteringMaterial(mu_s[1], mu_a[1], g, n[1]), label="middleLayer") + backLayer = Cuboid(w, w, t[2], material=ScatteringMaterial(mu_s[2], mu_a[2], g, n[2]), label="backLayer") layerStack = backLayer.stack(middleLayer, "front").stack(frontLayer, "front") layerStack.translateTo(Vector(0, 0, sum(t) / 2)) diff --git a/pytissueoptics/rayscattering/tests/display/testProfileFactory.py b/pytissueoptics/rayscattering/tests/display/testProfileFactory.py index 16410809..a49afb27 100644 --- a/pytissueoptics/rayscattering/tests/display/testProfileFactory.py +++ b/pytissueoptics/rayscattering/tests/display/testProfileFactory.py @@ -127,7 +127,7 @@ def testWhenCreateProfile_shouldSetProperName(self): profile = self.profileFactory.create( Direction.Z_POS, solidLabel="cube", surfaceLabel="top", surfaceEnergyLeaving=True ) - self.assertEqual("Energy profile along z of cube surface top (leaving)", profile.name) + self.assertEqual("DEPOSITION profile along z of cube surface top (leaving)", profile.name) def testGivenEmptyLogger_whenCreateProfile_shouldReturnEmptyProfile(self): emptyLogger = EnergyLogger(self.TEST_SCENE) diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py index b2edaeb4..a4e036c6 100644 --- a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py +++ b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py @@ -215,7 +215,7 @@ def testGivenALoggerPreviouslySaved_whenLoad_shouldLoadPreviousLoggerFromFile(se logger = EnergyLogger(self.TEST_SCENE) logger.load(filePath) - self.assertTrue(np.array_equal(previousLogger.getDataPoints(), logger.getDataPoints())) + self.assertTrue(np.array_equal(previousLogger.getRawDataPoints(), logger.getRawDataPoints())) self.assertEqual(previousLogger.info, logger.info) def testGivenALoggerPreviouslySaved_whenCreatingNewLoggerFromFile_shouldLoadPreviousLoggerFromFile(self): @@ -229,7 +229,7 @@ def testGivenALoggerPreviouslySaved_whenCreatingNewLoggerFromFile_shouldLoadPrev logger = EnergyLogger(self.TEST_SCENE, filePath) - self.assertTrue(np.array_equal(previousLogger.getDataPoints(), logger.getDataPoints())) + self.assertTrue(np.array_equal(previousLogger.getRawDataPoints(), logger.getRawDataPoints())) self.assertEqual(previousLogger.info, logger.info) def testGivenLoggerFromFile_shouldWarnIfLoadedWithDifferentScene(self): diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testPointCloudFactory.py b/pytissueoptics/rayscattering/tests/energyLogging/testPointCloudFactory.py index 37088ec9..9a107aa7 100644 --- a/pytissueoptics/rayscattering/tests/energyLogging/testPointCloudFactory.py +++ b/pytissueoptics/rayscattering/tests/energyLogging/testPointCloudFactory.py @@ -2,8 +2,9 @@ import numpy as np +from pytissueoptics import EnergyLogger, ScatteringScene from pytissueoptics.rayscattering.energyLogging import PointCloudFactory -from pytissueoptics.scene.logger import InteractionKey, Logger +from pytissueoptics.scene.logger import InteractionKey class TestPointCloudFactory(unittest.TestCase): @@ -55,7 +56,7 @@ def testWhenGetPointCloudOfSurfaces_shouldReturnPointCloudWithAllSurfacePoints(s self.assertIsNone(pointCloud.solidPoints) def testGivenEmptyLogger_whenGetPointCloud_shouldReturnEmptyPointCloud(self): - logger = Logger() + logger = EnergyLogger(scene=ScatteringScene([])) pointCloudFactory = PointCloudFactory(logger) pointCloud = pointCloudFactory.getPointCloud() @@ -63,7 +64,7 @@ def testGivenEmptyLogger_whenGetPointCloud_shouldReturnEmptyPointCloud(self): self.assertIsNone(pointCloud.surfacePoints) def _createTestLogger(self): - logger = Logger() + logger = EnergyLogger(scene=ScatteringScene([])) solidPointsA = np.array([[0.5, 1, 0, 0], [0.5, 1, 0, 0.1], [0.5, 1, 0, -0.1]]) solidPointsB = np.array([[0.5, -1, 0, 0], [0.5, -1, 0, 0.1], [0.5, -1, 0, -0.1]]) surfacePointsA = np.array([[1, 1, 0, 0.1], [-1, 1, 0, -0.1]]) diff --git a/pytissueoptics/rayscattering/tests/opencl/testCLPhotons.py b/pytissueoptics/rayscattering/tests/opencl/testCLPhotons.py index de0d8154..77b738d4 100644 --- a/pytissueoptics/rayscattering/tests/opencl/testCLPhotons.py +++ b/pytissueoptics/rayscattering/tests/opencl/testCLPhotons.py @@ -35,7 +35,7 @@ def testWhenPropagate_shouldPropagateUntilAllPhotonsHaveNoMoreEnergy(self): photons.propagate(IPP=IPP, verbose=False) - dataPoints = logger.getDataPoints() + dataPoints = logger.getRawDataPoints() totalWeightScattered = float(np.sum(dataPoints[:, 0])) # Roulette effect will result in total weight slightly different from N. self.assertAlmostEqual(N, totalWeightScattered, places=1) @@ -58,16 +58,16 @@ def testWhenPropagateInSolids_shouldLogEnergyWithCorrectInteractionKeys(self): photons.propagate(IPP=IPP, verbose=False) - frontSurfacePoints = logger.getDataPoints(InteractionKey("cube", "cube_front")) + frontSurfacePoints = logger.getRawDataPoints(InteractionKey("cube", "cube_front")) energyInput = -np.sum(frontSurfacePoints[:, 0]) # should be around 97% of total energy because of reflections - cubePoints = logger.getDataPoints(InteractionKey("cube")) + cubePoints = logger.getRawDataPoints(InteractionKey("cube")) energyScattered = np.sum(cubePoints[:, 0]) energyLeaving = 0 for surfaceLabel in logger.getStoredSurfaceLabels("cube"): if "front" in surfaceLabel: continue - surfacePoints = logger.getDataPoints(InteractionKey("cube", surfaceLabel)) + surfacePoints = logger.getRawDataPoints(InteractionKey("cube", surfaceLabel)) energyLeaving += np.sum(surfacePoints[:, 0]) self.assertAlmostEqual(energyInput, energyScattered + energyLeaving, places=2) @@ -88,6 +88,6 @@ def testWhenPropagateOnly1Photon_shouldPropagate(self): photons.propagate(IPP=IPP, verbose=False) - dataPoints = logger.getDataPoints() + dataPoints = logger.getRawDataPoints() totalWeightScattered = float(np.sum(dataPoints[:, 0])) self.assertAlmostEqual(N, totalWeightScattered, places=2) diff --git a/pytissueoptics/scene/logger/logger.py b/pytissueoptics/scene/logger/logger.py index 15431d3d..2e813795 100644 --- a/pytissueoptics/scene/logger/logger.py +++ b/pytissueoptics/scene/logger/logger.py @@ -3,7 +3,7 @@ import warnings from dataclasses import dataclass from enum import Enum -from typing import Dict, List, Optional, Union +from typing import Callable, Dict, List, Optional, Union import numpy as np @@ -16,6 +16,17 @@ class InteractionKey: solidLabel: str surfaceLabel: Optional[str] = None + @property + def volumetric(self) -> bool: + return self.surfaceLabel is None + + +TransformFn = Callable[[InteractionKey, Optional[np.ndarray]], Optional[np.ndarray]] + + +def _noTransform(_: InteractionKey, data: Optional[np.ndarray]) -> Optional[np.ndarray]: + return data + @dataclass class InteractionData: @@ -74,17 +85,17 @@ def logSegment(self, start: Vector, end: Vector, key: InteractionKey = None): self._appendData([start.x, start.y, start.z, end.x, end.y, end.z], DataType.SEGMENT, key) def logPointArray(self, array: np.ndarray, key: InteractionKey = None): - """'array' must be of shape (n, 3) where second axis is (x, y, z)""" + """'array' must be of shape (n, 3) where the second axis is (x, y, z)""" assert array.shape[1] == 3 and array.ndim == 2, "Point array must be of shape (n, 3)" self._appendData(array, DataType.POINT, key) def logDataPointArray(self, array: np.ndarray, key: InteractionKey): - """'array' must be of shape (n, 4) where second axis is (value, x, y, z)""" + """'array' must be of shape (n, 4) where the second axis is (value, x, y, z)""" assert array.shape[1] == 4 and array.ndim == 2, "Data point array must be of shape (n, 4)" self._appendData(array, DataType.DATA_POINT, key) def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None): - """'array' must be of shape (n, 6) where second axis is (x1, y1, z1, x2, y2, z2)""" + """'array' must be of shape (n, 6) where the second axis is (x1, y1, z1, x2, y2, z2)""" assert array.shape[1] == 6 and array.ndim == 2, "Segment array must be of shape (n, 6)" self._appendData(array, DataType.SEGMENT, key) @@ -115,22 +126,26 @@ def _validateKey(self, key: InteractionKey): def getPoints(self, key: InteractionKey = None) -> np.ndarray: return self._getData(DataType.POINT, key) - def getDataPoints(self, key: InteractionKey = None) -> np.ndarray: + def getRawDataPoints(self, key: InteractionKey = None) -> np.ndarray: + """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) where + the second axis is (value, x, y, z).""" return self._getData(DataType.DATA_POINT, key) def getSegments(self, key: InteractionKey = None) -> np.ndarray: return self._getData(DataType.SEGMENT, key) - def _getData(self, dataType: DataType, key: InteractionKey = None) -> Optional[np.ndarray]: + def _getData( + self, dataType: DataType, key: InteractionKey = None, transform: TransformFn = _noTransform + ) -> Optional[np.ndarray]: if key and key.solidLabel: if not self._keyExists(key): return None - container = getattr(self._data[key], dataType.value) - return container.getData() + container: ListArrayContainer = getattr(self._data[key], dataType.value) + return transform(key, container.getData()) else: container = ListArrayContainer() for interactionData in self._data.values(): - points = getattr(interactionData, dataType.value) + points: Optional[ListArrayContainer] = getattr(interactionData, dataType.value) if points is None: continue container.extend(points) diff --git a/pytissueoptics/scene/scene/scene.py b/pytissueoptics/scene/scene/scene.py index 838088c4..7e7b4085 100644 --- a/pytissueoptics/scene/scene/scene.py +++ b/pytissueoptics/scene/scene/scene.py @@ -149,6 +149,14 @@ def getMaterials(self) -> list: materials.append(material) return list(materials) + def getMaterial(self, solidLabel: str): + solid = self.getSolid(solidLabel) + if solid.isStack(): + layerSurfaceLabel = [s for s in solid.getLayerSurfaceLabels(solidLabel) if INTERFACE_KEY not in s][0] + return solid.getEnvironment(layerSurfaceLabel).material + else: + return solid.getEnvironment().material + def getBoundingBox(self) -> Optional[BoundingBox]: if len(self._solids) == 0: return None diff --git a/pytissueoptics/scene/tests/logger/testLogger.py b/pytissueoptics/scene/tests/logger/testLogger.py index d4d127c4..79dd2e22 100644 --- a/pytissueoptics/scene/tests/logger/testLogger.py +++ b/pytissueoptics/scene/tests/logger/testLogger.py @@ -17,7 +17,7 @@ def testGivenNewLogger_shouldBeEmpty(self): logger = Logger() self.assertIsNone(logger.getPoints()) - self.assertIsNone(logger.getDataPoints()) + self.assertIsNone(logger.getRawDataPoints()) self.assertIsNone(logger.getSegments()) def testWhenLogNewPoint_shouldAddPointToTheLoggedPoints(self): @@ -45,15 +45,15 @@ def testWhenLogNewDataPoint_shouldAddDataPointToTheLoggedDataPoints(self): logger.logDataPoint(10, Vector(2, 0, 0), self.INTERACTION_KEY) - self.assertEqual(3, len(logger.getDataPoints())) - self.assertTrue(np.array_equal([10, 2, 0, 0], logger.getDataPoints()[-1])) + self.assertEqual(3, len(logger.getRawDataPoints())) + self.assertTrue(np.array_equal([10, 2, 0, 0], logger.getRawDataPoints()[-1])) def testWhenLogDataPointArray_shouldAddAllDataPointsToTheLoggedDataPoints(self): logger = Logger() logger.logDataPointArray(np.array([[2, 0, 0, 0], [1, 1, 0, 0]]), self.INTERACTION_KEY) - self.assertEqual(2, len(logger.getDataPoints())) - self.assertTrue(np.array_equal([1, 1, 0, 0], logger.getDataPoints()[-1])) + self.assertEqual(2, len(logger.getRawDataPoints())) + self.assertTrue(np.array_equal([1, 1, 0, 0], logger.getRawDataPoints()[-1])) def testWhenLogNewSegment_shouldAddSegmentToTheLoggedSegments(self): logger = Logger() @@ -127,7 +127,7 @@ def testWhenGetSpecificDataTypeWithEmptyKey_shouldReturnAllDataOfThisType(self): logger.logDataPoint(0, Vector(0, 0, 0), self.INTERACTION_KEY) logger.logPoint(Vector(0, 0, 0), InteractionKey("another solid", "another surface")) - self.assertEqual(1, len(logger.getDataPoints())) + self.assertEqual(1, len(logger.getRawDataPoints())) def testWhenGetDataWithEmptyKey_shouldReturnAllData(self): anotherKey = InteractionKey(self.SOLID_LABEL, "another surface") @@ -157,7 +157,7 @@ def testGivenEmptyLogger_whenGetData_shouldReturnNone(self): self.assertIsNone(logger.getPoints()) self.assertIsNone(logger.getSegments()) - self.assertIsNone(logger.getDataPoints()) + self.assertIsNone(logger.getRawDataPoints()) def testWhenGetSolidLabels_shouldReturnAListOfUniqueSolidLabels(self): solidLabel1 = "A label" diff --git a/pytissueoptics/scene/viewer/mayavi/mayavi3DViewer.py b/pytissueoptics/scene/viewer/mayavi/mayavi3DViewer.py index 22da6b02..fbe34bea 100644 --- a/pytissueoptics/scene/viewer/mayavi/mayavi3DViewer.py +++ b/pytissueoptics/scene/viewer/mayavi/mayavi3DViewer.py @@ -67,7 +67,7 @@ def addLogger( ): self.addPoints(logger.getPoints(), colormap=colormap, reverseColormap=reverseColormap, scale=pointScale) self.addDataPoints( - logger.getDataPoints(), + logger.getRawDataPoints(), colormap=colormap, reverseColormap=reverseColormap, scale=dataPointScale,