Merge pull request #138 from DCC-Lab/jlb/fluence

Add fluence rate option.

Author: JLBegin

pytissueoptics/examples/benchmarks/cube60.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "MCX Homogeneous cube"

pytissueoptics/examples/benchmarks/cubesph60b.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "MCX Sphere"

pytissueoptics/examples/benchmarks/skinvessel.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "Skin vessel"

pytissueoptics/examples/benchmarks/sphshells.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "MCX Spherical shells"

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)

pytissueoptics/examples/rayscattering/ex02.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "Propagation in an Infinite Medium"

pytissueoptics/examples/rayscattering/ex03.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = (

pytissueoptics/examples/rayscattering/ex04.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "Custom layer stack"

pytissueoptics/examples/rayscattering/ex05.py

@@ -1,4 +1,5 @@
 import env  # noqa: F401
+
 from pytissueoptics import *  # noqa: F403
 
 TITLE = "Sphere inside a cube"

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",

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()

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:

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,16 +203,21 @@ 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):
         stats = Stats(self._logger)
         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)