Scene Module

Cuboid

Bases: Solid

Also known as the Right Rectangular Prism, the Cuboid is defined by its width (a, b, c) in each axis (x, y, z) respectively.

The position refers to the vector from global origin to its centroid. The generated mesh will be divided into the following subgroups: Left (-x), Right (+x), Bottom (-y), Top (+y), Front (-z), Back (+z).

Source code in pytissueoptics/scene/solids/cuboid.py

class Cuboid(Solid):
    """
    Also known as the Right Rectangular Prism, the Cuboid is defined by its
    width (a, b, c) in each axis (x, y, z) respectively.

    The position refers to the vector from global origin to its centroid.
    The generated mesh will be divided into the following subgroups:
    Left (-x), Right (+x), Bottom (-y), Top (+y), Front (-z), Back (+z).
    """

    def __init__(
        self,
        a: float,
        b: float,
        c: float,
        vertices: List[Vertex] = None,
        position: Vector = Vector(0, 0, 0),
        surfaces: SurfaceCollection = None,
        material=None,
        label: str = "cuboid",
        primitive: str = primitives.DEFAULT,
        labelOverride=True,
    ):
        self.shape = [a, b, c]

        if not vertices:
            vertices = [
                Vertex(-a / 2, -b / 2, c / 2),
                Vertex(a / 2, -b / 2, c / 2),
                Vertex(a / 2, b / 2, c / 2),
                Vertex(-a / 2, b / 2, c / 2),
                Vertex(-a / 2, -b / 2, -c / 2),
                Vertex(a / 2, -b / 2, -c / 2),
                Vertex(a / 2, b / 2, -c / 2),
                Vertex(-a / 2, b / 2, -c / 2),
            ]

        super().__init__(vertices, position, surfaces, material, label, primitive, labelOverride=labelOverride)

    def _computeTriangleMesh(self):
        V = self._vertices
        self._surfaces.add("left", [Triangle(V[4], V[0], V[3]), Triangle(V[3], V[7], V[4])])
        self._surfaces.add("right", [Triangle(V[1], V[5], V[6]), Triangle(V[6], V[2], V[1])])
        self._surfaces.add("bottom", [Triangle(V[4], V[5], V[1]), Triangle(V[1], V[0], V[4])])
        self._surfaces.add("top", [Triangle(V[3], V[2], V[6]), Triangle(V[6], V[7], V[3])])
        self._surfaces.add("front", [Triangle(V[5], V[4], V[7]), Triangle(V[7], V[6], V[5])])
        self._surfaces.add("back", [Triangle(V[0], V[1], V[2]), Triangle(V[2], V[3], V[0])])

    def _computeQuadMesh(self):
        V = self._vertices
        self._surfaces.add("left", [Quad(V[4], V[0], V[3], V[7])])
        self._surfaces.add("right", [Quad(V[1], V[5], V[6], V[2])])
        self._surfaces.add("bottom", [Quad(V[4], V[5], V[1], V[0])])
        self._surfaces.add("top", [Quad(V[3], V[2], V[6], V[7])])
        self._surfaces.add("front", [Quad(V[5], V[4], V[7], V[6])])
        self._surfaces.add("back", [Quad(V[0], V[1], V[2], V[3])])

    def stack(self, other: "Cuboid", onSurface: str = "top", stackLabel="CuboidStack") -> "Cuboid":
        """
        Basic implementation for stacking cuboids along an axis.

        For example, stacking on 'top' will move the other cuboid on top of this cuboid. They will now share
         the same mesh at the interface and inside/outside materials at the interface will be properly defined.
         This will return a new cuboid that represents the stack, with a new 'interface<i>' surface group.

        Limitations:
            - Requires cuboids with the same shape except along the stack axis.
            - Cannot stack another stack unless it is along its stacked axis (ill-defined interface material).
            - Expected behavior not guaranteed for pre-rotated cuboids.
            - Stacked cuboids will lose reference to their initial stack surface (not a closed solid anymore).
                Use the returned cuboid stack.
        """
        stacker = CuboidStacker()
        stackResult = stacker.stack(onCuboid=self, otherCuboid=other, onSurface=onSurface)
        return Cuboid._fromStackResult(stackResult, label=stackLabel)

    @classmethod
    def _fromStackResult(cls, stackResult: StackResult, label: str) -> "Cuboid":
        # subtracting stackCentroid from all vertices because solid creation will translate back to position.
        for vertex in stackResult.vertices:
            vertex.subtract(stackResult.position)

        cuboid = Cuboid(
            *stackResult.shape,
            position=stackResult.position,
            vertices=stackResult.vertices,
            surfaces=stackResult.surfaces,
            label=label,
            primitive=stackResult.primitive,
            labelOverride=False,
        )
        cuboid._layerLabels = stackResult.layerLabels
        return cuboid

    def contains(self, *vertices: Vector) -> bool:
        relativeVertices = [vertex - self.position for vertex in vertices]
        relativeVertices = self._applyInverseRotation(relativeVertices)

        relativeVertices = np.asarray([vertex.array for vertex in relativeVertices])
        bounds = [s / 2 for s in self.shape]
        if np.any(np.abs(relativeVertices) >= bounds):
            return False

        if len(vertices) == 1:
            return True

        if self.isStack():
            # At this point, all vertices are contained in the bounding box of the stack.
            # We just need to make sure they are contained in a single layer.
            return self._aSingleLayerContains(relativeVertices)
        return True

    def _aSingleLayerContains(self, relativeVertices: np.ndarray) -> bool:
        for layerLabel in self._layerLabels:
            bbox = self._getLayerBBox(layerLabel)
            bboxShape = [bbox.xWidth, bbox.yWidth, bbox.zWidth]
            bounds = [s / 2 for s in bboxShape]
            layerOffset = bbox.center - self.position
            layerRelativeVertices = relativeVertices - np.asarray(layerOffset.array)
            if np.all(np.abs(layerRelativeVertices) < bounds):
                return True
        return False

    def _getLayerBBox(self, layerLabel: str) -> BoundingBox:
        polygons = self._getLayerPolygons(layerLabel)
        bbox = BoundingBox.fromPolygons(polygons)
        return bbox

    def _getLayerPolygons(self, layerLabel: str) -> List[Polygon]:
        layerSurfaceLabels = self._layerLabels[layerLabel]
        polygons = []
        for surfaceLabel in layerSurfaceLabels:
            polygons.extend(self._surfaces.getPolygons(surfaceLabel))
        return polygons

    def _geometryParams(self) -> dict:
        return {"shape": self.shape}

    def __hash__(self):
        baseHash = super().__hash__()
        if not self.isStack():
            return baseHash

        materials = set()
        for surfaceLabel in self.surfaceLabels:
            material = self.getPolygons(surfaceLabel)[0].insideEnvironment.material
            materials.add(material)
        materialsHash = hash(tuple(materials))
        return hash((baseHash, materialsHash))

stack(other, onSurface='top', stackLabel='CuboidStack')

Basic implementation for stacking cuboids along an axis.

For example, stacking on 'top' will move the other cuboid on top of this cuboid. They will now share the same mesh at the interface and inside/outside materials at the interface will be properly defined. This will return a new cuboid that represents the stack, with a new 'interface' surface group.

Limitations

• Requires cuboids with the same shape except along the stack axis.
• Cannot stack another stack unless it is along its stacked axis (ill-defined interface material).
• Expected behavior not guaranteed for pre-rotated cuboids.
• Stacked cuboids will lose reference to their initial stack surface (not a closed solid anymore). Use the returned cuboid stack.

Source code in pytissueoptics/scene/solids/cuboid.py

def stack(self, other: "Cuboid", onSurface: str = "top", stackLabel="CuboidStack") -> "Cuboid":
    """
    Basic implementation for stacking cuboids along an axis.

    For example, stacking on 'top' will move the other cuboid on top of this cuboid. They will now share
     the same mesh at the interface and inside/outside materials at the interface will be properly defined.
     This will return a new cuboid that represents the stack, with a new 'interface<i>' surface group.

    Limitations:
        - Requires cuboids with the same shape except along the stack axis.
        - Cannot stack another stack unless it is along its stacked axis (ill-defined interface material).
        - Expected behavior not guaranteed for pre-rotated cuboids.
        - Stacked cuboids will lose reference to their initial stack surface (not a closed solid anymore).
            Use the returned cuboid stack.
    """
    stacker = CuboidStacker()
    stackResult = stacker.stack(onCuboid=self, otherCuboid=other, onSurface=onSurface)
    return Cuboid._fromStackResult(stackResult, label=stackLabel)

Cylinder

Bases: Solid

Source code in pytissueoptics/scene/solids/cylinder.py

class Cylinder(Solid):
    def __init__(
        self,
        radius: float = 1,
        length: float = 1,
        u: int = 32,
        v: int = 3,
        s: int = 2,
        position: Vector = Vector(0, 0, 0),
        material=None,
        primitive: str = primitives.DEFAULT,
        label: str = "cylinder",
        smooth=True,
    ):
        """
        Default cylinder orientation will be along the z axis. The front face towards the negative z axis and the back
        face towards the positive z axis. Position refers to its centroid. Available surfaces are "front", "lateral" and
        "back".
        """
        self._radius = radius
        self._length = length
        if u < 3:
            raise ValueError("u must be >= 3")
        if length != 0 and v < 1:
            raise ValueError("v must be >= 1 for non-zero length")
        assert radius > 0 and length >= 0
        self._u = u
        self._v = v
        self._s = s
        self._frontCenter = Vertex(0, 0, 0)
        self._backCenter = Vertex(0, 0, length)
        self._minRadius = math.cos(math.pi / self._u) * self._radius
        self._angularStep = 2 * math.pi / self._u
        self._lateralStep = self._length / self._v
        self._radialStep = self._radius / self._s

        super().__init__(
            position=position, material=material, primitive=primitive, vertices=[], smooth=smooth, label=label
        )
        self.translateBy(Vector(0, 0, -length / 2))
        self._position += Vector(0, 0, length / 2)

    @property
    def direction(self) -> Vector:
        direction = self._backCenter - self._frontCenter
        direction.normalize()
        return direction

    def _computeTriangleMesh(self):
        frontLayers, lateralLayers, backLayers = self._computeVerticesOfLayers()

        frontLayers.insert(0, [self._frontCenter])
        frontLayers.append(lateralLayers[0])
        backLayers.insert(0, lateralLayers[-1])
        backLayers.append([self._backCenter])
        self._vertices.extend([self._frontCenter, self._backCenter])

        self._surfaces.add("front", self._getSurfaceTriangles(frontLayers))
        self._surfaces.add("lateral", self._getSurfaceTriangles(lateralLayers))
        self._surfaces.add("back", self._getSurfaceTriangles(backLayers))

    def _computeVerticesOfLayers(self) -> tuple:
        v = self._v if self._length != 0 else 0
        frontLayers = self._computeSectionVertices(lateralSteps=[0], radialSteps=list(range(1, self._s)))
        lateralLayers = self._computeSectionVertices(lateralSteps=list(range(v + 1)), radialSteps=[self._s])
        backLayers = self._computeSectionVertices(lateralSteps=[v], radialSteps=list(range(self._s - 1, 0, -1)))
        return frontLayers, lateralLayers, backLayers

    def _computeSectionVertices(self, lateralSteps: List[int], radialSteps: List[int]):
        verticesLayers = []
        for k in radialSteps:
            for j in lateralSteps:
                layer = []
                for i in range(self._u):
                    layer.append(self._createVertex(i, j, k))
                verticesLayers.append(layer)
                self._vertices.extend(layer)
        return verticesLayers

    def _createVertex(self, i: int, j: int, k: int) -> Vertex:
        shrinkFactor = self._getShrinkFactor(self._lateralStep * j)
        radiusFactor = k * self._radialStep / self._radius
        if shrinkFactor != 1:
            pass  # This is for cones only. In this case we want to keep linear sampling of radius values.
        else:
            # For lenses, we usually want a uniform mesh after the curve transform, but this transform tends to push
            # vertices outwards (particularly at low radius). To prevent a low mesh resolution in the center, we need
            # to increase the sampling around the center beforehand by forcing smaller radiusFactor values.
            radiusFactor = radiusFactor**2
        r = self._radius * radiusFactor * shrinkFactor
        x = r * math.cos(i * self._angularStep)
        y = r * math.sin(i * self._angularStep)
        z = j * self._lateralStep
        return Vertex(x, y, z)

    def _getSurfaceTriangles(self, verticesLayers: List[List[Vertex]]) -> List[Triangle]:
        triangles = []
        for i in range(len(verticesLayers) - 1):
            currentGroup = verticesLayers[i]
            nextGroup = verticesLayers[i + 1]

            if len(currentGroup) == 1:
                triangles.extend(self._getPeakTriangles(self._frontCenter, nextGroup[::-1]))
                continue
            if len(nextGroup) == 1:
                triangles.extend(self._getPeakTriangles(self._backCenter, currentGroup))
                continue

            for j in range(self._u):
                nextIndex = (j + 1) % self._u
                triangles.append(Triangle(currentGroup[j], nextGroup[nextIndex], nextGroup[j]))
                triangles.append(Triangle(currentGroup[j], currentGroup[nextIndex], nextGroup[nextIndex]))
        return triangles

    def _getPeakTriangles(self, peakVertex: Vertex, ringVertices: List[Vertex]) -> List[Triangle]:
        triangles = []
        for j in range(self._u):
            nextIndex = (j + 1) % self._u
            triangles.append(Triangle(peakVertex, ringVertices[j], ringVertices[nextIndex]))
        return triangles

    def _computeQuadMesh(self):
        raise NotImplementedError("Quad mesh not implemented for Cylinder")

    def contains(self, *vertices: Vector) -> bool:
        direction = self.direction
        basePosition = self._position - direction * self._length / 2
        for vertex in vertices:
            localPoint = vertex - basePosition
            alongCylinder = direction.dot(localPoint)
            if alongCylinder < 0 or alongCylinder > self._length:
                return False
            radiusCheck = self._minRadiusAtHeightAlong(alongCylinder)
            radialComponent = localPoint - direction * alongCylinder
            radialDistanceFromBase = radialComponent.getNorm()
            if radialDistanceFromBase > radiusCheck:
                return False
        return True

    def _minRadiusAtHeightAlong(self, heightAlong: float) -> float:
        shrinkFactor = self._getShrinkFactor(heightAlong)
        return self._minRadius * shrinkFactor

    @staticmethod
    def _getShrinkFactor(heightAlong: float) -> float:
        return 1

    def smooth(self, surfaceLabel: str = None, reset: bool = True):
        if self._u < 16:
            warnings.warn("Smoothing a cylinder with less than 16 sides (u < 16) may result in intersection errors.")
        if surfaceLabel:
            return super(Cylinder, self).smooth(surfaceLabel, reset)
        self.smooth("lateral", reset=True)

    def __hash__(self):
        materialHash = hash(self._material) if self._material else 0
        propertyHash = hash((self._radius, self._length, self._frontCenter, self._backCenter))
        return hash((materialHash, propertyHash))

    def _geometryParams(self) -> dict:
        return {
            "radius": self._radius,
            "length": self._length,
            "u": self._u,
            "v": self._v,
            "s": self._s,
        }

__init__(radius=1, length=1, u=32, v=3, s=2, position=Vector(0, 0, 0), material=None, primitive=primitives.DEFAULT, label='cylinder', smooth=True)

Default cylinder orientation will be along the z axis. The front face towards the negative z axis and the back face towards the positive z axis. Position refers to its centroid. Available surfaces are "front", "lateral" and "back".

Source code in pytissueoptics/scene/solids/cylinder.py

def __init__(
    self,
    radius: float = 1,
    length: float = 1,
    u: int = 32,
    v: int = 3,
    s: int = 2,
    position: Vector = Vector(0, 0, 0),
    material=None,
    primitive: str = primitives.DEFAULT,
    label: str = "cylinder",
    smooth=True,
):
    """
    Default cylinder orientation will be along the z axis. The front face towards the negative z axis and the back
    face towards the positive z axis. Position refers to its centroid. Available surfaces are "front", "lateral" and
    "back".
    """
    self._radius = radius
    self._length = length
    if u < 3:
        raise ValueError("u must be >= 3")
    if length != 0 and v < 1:
        raise ValueError("v must be >= 1 for non-zero length")
    assert radius > 0 and length >= 0
    self._u = u
    self._v = v
    self._s = s
    self._frontCenter = Vertex(0, 0, 0)
    self._backCenter = Vertex(0, 0, length)
    self._minRadius = math.cos(math.pi / self._u) * self._radius
    self._angularStep = 2 * math.pi / self._u
    self._lateralStep = self._length / self._v
    self._radialStep = self._radius / self._s

    super().__init__(
        position=position, material=material, primitive=primitive, vertices=[], smooth=smooth, label=label
    )
    self.translateBy(Vector(0, 0, -length / 2))
    self._position += Vector(0, 0, length / 2)

Ellipsoid

Bases: Solid

We take the unit sphere, then calculate the theta, phi position of each vertex (with ISO mathematical convention). Then we apply the ellipsoid formula in the spherical coordinate to isolate the component R. We then calculate the difference the ellipsoid would with the unit sphere for this theta,phi and then .add() or .subtract() the corresponding vector.

Source code in pytissueoptics/scene/solids/ellipsoid.py

class Ellipsoid(Solid):
    """
    We take the unit sphere, then calculate the theta, phi position of each vertex (with ISO mathematical
    convention). Then we apply the ellipsoid formula in the spherical coordinate to isolate the component R.
    We then calculate the difference the ellipsoid would with the unit sphere for this theta,phi and
    then .add() or .subtract() the corresponding vector.
    """

    def __init__(
        self,
        a: float = 1,
        b: float = 1,
        c: float = 1,
        order: int = 3,
        position: Vector = Vector(0, 0, 0),
        material=None,
        label: str = "ellipsoid",
        primitive: str = primitives.DEFAULT,
        smooth: bool = True,
    ):
        self._a = a
        self._b = b
        self._c = c
        self._order = order

        super().__init__(
            position=position, material=material, label=label, primitive=primitive, vertices=[], smooth=smooth
        )

    def _computeTriangleMesh(self):
        """
        The most precise sphere approximation is the IcoSphere, which is generated from the platonic solid,
        the Icosahedron. It is built with 20 equilateral triangles with exactly the same angle between each.
        From Euler's method to generate the vertex for the icosahedron, we cross 3 perpendicular planes,
        with length=2 and width=2*phi. Joining adjacent vertices will produce the Icosahedron.

        From the Icosahedron, we can split each face in 4 triangles, in a recursive manner, to obtain an IcoSphere.
        The method goes as follows:
        1 - Find each mid-point between two connecting vertices on a triangle
        2 - Normalize those new points to project them onto the unit sphere.
        3 - Connect the new vertices in a way to make 4 new triangles.
        4 - Do these steps for each triangle (This will lead to redundant calculation, I am aware)
        5 - Replace the old surfaces by the new surfaces
        """

        self._computeFirstOrderTriangleMesh()

        for i in range(0, self._order):
            self._computeNextOrderTriangleMesh()

        self._setVerticesPositionsFromCenter()

    def _computeFirstOrderTriangleMesh(self):
        phi = (1.0 + 5.0 ** (1 / 2)) / 2.0
        xyPlaneVertices = [Vertex(-1, phi, 0), Vertex(1, phi, 0), Vertex(-1, -phi, 0), Vertex(1, -phi, 0)]
        yzPlaneVertices = [Vertex(0, -1, phi), Vertex(0, 1, phi), Vertex(0, -1, -phi), Vertex(0, 1, -phi)]
        xzPlaneVertices = [Vertex(phi, 0, -1), Vertex(phi, 0, 1), Vertex(-phi, 0, -1), Vertex(-phi, 0, 1)]
        self._vertices = [*xyPlaneVertices, *yzPlaneVertices, *xzPlaneVertices]
        V = self._vertices

        self._surfaces.add(
            "ellipsoid",
            [
                Triangle(V[0], V[11], V[5]),
                Triangle(V[0], V[5], V[1]),
                Triangle(V[0], V[1], V[7]),
                Triangle(V[0], V[7], V[10]),
                Triangle(V[0], V[10], V[11]),
                Triangle(V[1], V[5], V[9]),
                Triangle(V[5], V[11], V[4]),
                Triangle(V[11], V[10], V[2]),
                Triangle(V[10], V[7], V[6]),
                Triangle(V[7], V[1], V[8]),
                Triangle(V[3], V[9], V[4]),
                Triangle(V[3], V[4], V[2]),
                Triangle(V[3], V[2], V[6]),
                Triangle(V[3], V[6], V[8]),
                Triangle(V[3], V[8], V[9]),
                Triangle(V[4], V[9], V[5]),
                Triangle(V[2], V[4], V[11]),
                Triangle(V[6], V[2], V[10]),
                Triangle(V[8], V[6], V[7]),
                Triangle(V[9], V[8], V[1]),
            ],
        )

    def _computeNextOrderTriangleMesh(self):
        newPolygons = []
        self._verticesCache = {}
        for j, polygon in enumerate(self.getPolygons()):
            ai = self._createMidVertex(polygon.vertices[0], polygon.vertices[1])
            bi = self._createMidVertex(polygon.vertices[1], polygon.vertices[2])
            ci = self._createMidVertex(polygon.vertices[2], polygon.vertices[0])
            newVertices = [ai, bi, ci]
            for i, vertex in enumerate(newVertices):
                vHash = hash2((vertex.x, vertex.y, vertex.z))
                if vHash in self._verticesCache:
                    newVertices[i] = self._verticesCache[vHash]
                    self._verticesCache.pop(vHash)
                else:
                    self._vertices.append(vertex)
                    self._verticesCache[vHash] = vertex

            newPolygons.append(Triangle(polygon.vertices[0], newVertices[0], newVertices[2]))
            newPolygons.append(Triangle(polygon.vertices[1], newVertices[1], newVertices[0]))
            newPolygons.append(Triangle(polygon.vertices[2], newVertices[2], newVertices[1]))
            newPolygons.append(Triangle(newVertices[0], newVertices[1], newVertices[2]))

        self._surfaces.setPolygons("ellipsoid", newPolygons)

    @staticmethod
    def _createMidVertex(p1, p2):
        middle = Vertex((p1.x + p2.x) / 2, (p1.y + p2.y) / 2, (p1.z + p2.z) / 2)
        return middle

    def _setVerticesPositionsFromCenter(self):
        """
        The Ellipsoid parametric equation goes as: x^2/a^2 + y^2/b^2 + z^2/c^2 =1
        A Sphere is just an ellipsoid with a = b = c.
        Bringing (x, y, z) -> (theta, phi, r) we can simply take the unit sphere and stretch it,
        since the equation becomes as follow:

        r^2.cos^2(theta).sin^2(phi)/a^2 + r^2.sin^2(theta).sin^2(phi)/b^2  + r^2.cos^2(phi)/c^2 = 1
        """
        for vertex in self._vertices:
            vertex.normalize()
            r = self._radiusTowards(vertex)
            distanceFromUnitSphere = r - 1.0
            vertex.add(vertex * distanceFromUnitSphere)
        self.surfaces.resetNormals()

    @staticmethod
    def _findThetaPhi(vertex: Vertex):
        phi = math.acos(vertex.z / (vertex.x**2 + vertex.y**2 + vertex.z**2))
        theta = 0
        if vertex.x == 0.0:
            if vertex.y > 0.0:
                theta = math.pi / 2

            elif vertex.y < 0.0:
                theta = -math.pi / 2

        elif vertex.x > 0.0:
            theta = math.atan(vertex.y / vertex.x)

        elif vertex.x < 0.0:
            if vertex.y >= 0.0:
                theta = math.atan(vertex.y / vertex.x) + math.pi

            elif vertex.y < 0.0:
                theta = math.atan(vertex.y / vertex.x) - math.pi

        return theta, phi

    def _radiusTowards(self, vertex):
        theta, phi = self._findThetaPhi(vertex)
        return math.sqrt(
            1
            / (
                (math.cos(theta) ** 2 * math.sin(phi) ** 2) / self._a**2
                + (math.sin(theta) ** 2 * math.sin(phi) ** 2) / self._b**2
                + math.cos(phi) ** 2 / self._c**2
            )
        )

    def _computeQuadMesh(self):
        raise NotImplementedError

    def contains(self, *vertices: Vector) -> bool:
        """Only returns true if all vertices are inside the minimum radius of the ellipsoid
        towards each vertex direction (more restrictive with low order ellipsoids)."""
        relativeVertices = [vertex - self.position for vertex in vertices]
        relativeVertices = self._applyInverseRotation(relativeVertices)
        relativeVerticesArray = np.asarray([vertex.array for vertex in relativeVertices])

        for relativeVertexArray in relativeVerticesArray:
            relativeVertex = Vertex(*relativeVertexArray)
            vertexRadius = relativeVertex.getNorm()
            relativeVertex.normalize()
            if vertexRadius == 0:
                continue
            minRadius = self._getMinimumRadiusTowards(relativeVertex)
            if vertexRadius >= minRadius:
                return False
        return True

    def _getRadiusError(self) -> float:
        aPolygon = self.surfaces.getPolygons()[0]
        centerVertex = Vertex(0, 0, 0)
        for vertex in aPolygon.vertices:
            centerVertex.add(vertex)
        centerVertex.divide(len(aPolygon.vertices))
        centerVertex.subtract(self.position)

        localMinimumRadius = centerVertex.getNorm()
        localTrueRadius = self._radiusTowards(centerVertex)
        return abs(localTrueRadius - localMinimumRadius) / localTrueRadius

    def _getMinimumRadiusTowards(self, vertex) -> float:
        return (1 - self._getRadiusError()) * self._radiusTowards(vertex)

    def _geometryParams(self) -> dict:
        return {
            "radius_a": self._a,
            "radius_b": self._b,
            "radius_c": self._c,
            "order": self._order,
        }

contains(*vertices)

Only returns true if all vertices are inside the minimum radius of the ellipsoid towards each vertex direction (more restrictive with low order ellipsoids).

Source code in pytissueoptics/scene/solids/ellipsoid.py

def contains(self, *vertices: Vector) -> bool:
    """Only returns true if all vertices are inside the minimum radius of the ellipsoid
    towards each vertex direction (more restrictive with low order ellipsoids)."""
    relativeVertices = [vertex - self.position for vertex in vertices]
    relativeVertices = self._applyInverseRotation(relativeVertices)
    relativeVerticesArray = np.asarray([vertex.array for vertex in relativeVertices])

    for relativeVertexArray in relativeVerticesArray:
        relativeVertex = Vertex(*relativeVertexArray)
        vertexRadius = relativeVertex.getNorm()
        relativeVertex.normalize()
        if vertexRadius == 0:
            continue
        minRadius = self._getMinimumRadiusTowards(relativeVertex)
        if vertexRadius >= minRadius:
            return False
    return True

Loader

Base class to manage the conversion between files and Scene() or Solid() from various types of files.

Source code in pytissueoptics/scene/loader/loader.py

class Loader:
    """
    Base class to manage the conversion between files and Scene() or Solid() from
    various types of files.
    """

    def __init__(self):
        self._filepath: str = ""
        self._fileExtension: str = ""
        self._parser = None

    def load(self, filepath: str, showProgress: bool = True) -> List[Solid]:
        self._filepath = filepath
        self._fileExtension = self._getFileExtension()
        self._selectParser(showProgress)
        return self._convert(showProgress)

    def _getFileExtension(self) -> str:
        return pathlib.Path(self._filepath).suffix

    def _selectParser(self, showProgress: bool = True):
        ext = self._fileExtension
        if ext == ".obj":
            self._parser = OBJParser(self._filepath, showProgress)
        else:
            raise NotImplementedError("This format is not supported.")

    def _convert(self, showProgress: bool = True) -> List[Solid]:
        vertices = []
        for vertex in self._parser.vertices:
            vertices.append(Vertex(*vertex))

        totalProgressBarLength = 0
        for objectName, _object in self._parser.objects.items():
            totalProgressBarLength += len(_object.surfaces.items())
        pbar = progressBar(
            total=totalProgressBarLength,
            desc="Converting File '{}'".format(self._filepath.split("/")[-1]),
            unit="surfaces",
            disable=not showProgress,
        )

        solids = []
        for objectName, _object in self._parser.objects.items():
            surfaces = SurfaceCollection()
            for surfaceLabel, surface in _object.surfaces.items():
                surfaces.add(surfaceLabel, self._convertSurfaceToTriangles(surface, vertices))
                pbar.update(1)
            solids.append(
                Solid(
                    position=Vector(0, 0, 0),
                    vertices=vertices,
                    surfaces=surfaces,
                    primitive=primitives.POLYGON,
                    label=objectName,
                )
            )

        pbar.close()
        return solids

    @staticmethod
    def _convertSurfaceToTriangles(surface: ParsedSurface, vertices: List[Vertex]) -> List[Triangle]:
        """Converting to triangles only since loaded polygons are often not planar."""
        triangles = []
        for polygonIndices in surface.polygons:
            polygonVertices = [vertices[i] for i in polygonIndices]
            for i in range(len(polygonVertices) - 2):
                triangles.append(Triangle(polygonVertices[0], polygonVertices[i + 1], polygonVertices[i + 2]))
        return triangles

Logger

Source code in pytissueoptics/scene/logger/logger.py

class Logger:
    DEFAULT_LOGGER_PATH = "simulation.log"

    def __init__(self, fromFilepath: str = None):
        self._data: Dict[InteractionKey, InteractionData] = {}
        self.info: dict = {}
        self._filepath = None
        self._labels = {}

        if fromFilepath:
            self.load(fromFilepath)

    def getSeenSolidLabels(self) -> List[str]:
        """Returns a list of all solid labels that have been logged in the past
        even if the data was discarded."""
        return list(self._labels.keys())

    def getSeenSurfaceLabels(self, solidLabel: str) -> List[str]:
        """Returns a list of all surface labels that have been logged in the past
        for the given solid even if the data was discarded."""
        return self._labels[solidLabel]

    def getStoredSolidLabels(self) -> List[str]:
        """Returns a list of all solid labels that are currently stored in the logger."""
        return list(set(key.solidLabel for key in self._data.keys()))

    def getStoredSurfaceLabels(self, solidLabel: str) -> List[str]:
        """Returns a list of all surface labels that are currently stored in the logger."""
        return [
            key.surfaceLabel
            for key in self._data.keys()
            if key.solidLabel == solidLabel and key.surfaceLabel is not None
        ]

    def logPoint(self, point: Vector, key: InteractionKey = None):
        self._appendData([point.x, point.y, point.z], DataType.POINT, key)

    def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None):
        dataPoint = [value, *position.array]
        if ID is not None:
            dataPoint.append(ID)
        self._appendData(dataPoint, DataType.DATA_POINT, key)

    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 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) or (n, 5) where the second axis is (value, x, y, z)
        or (value, x, y, z, photonID). The photonID column is optional."""
        assert array.shape[1] in [4, 5] and array.ndim == 2, "Data point array must be of shape (n, 4) or (n, 5)"
        self._appendData(array, DataType.DATA_POINT, key)

    def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None):
        """'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)

    def _appendData(self, data: Union[List, np.ndarray], dataType: DataType, key: InteractionKey = None):
        if key is None:
            key = InteractionKey(None, None)
        self._validateKey(key)
        previousData = getattr(self._data[key], dataType.value)
        if previousData is None:
            previousData = ListArrayContainer()
            previousData.append(data)
            setattr(self._data[key], dataType.value, previousData)
        else:
            previousData.append(data)

    def _validateKey(self, key: InteractionKey):
        if key not in self._data:
            self._data[key] = InteractionData()
        if key.solidLabel is None:
            return
        if key.solidLabel not in self._labels:
            self._labels[key.solidLabel] = []
        if key.surfaceLabel is None:
            return
        if key.surfaceLabel not in self._labels[key.solidLabel]:
            self._labels[key.solidLabel].append(key.surfaceLabel)

    def getPoints(self, key: InteractionKey = None) -> np.ndarray:
        return self._getData(DataType.POINT, key)

    def getRawDataPoints(self, key: InteractionKey = None) -> np.ndarray:
        """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) or (n, 5) where
        the second axis is (value, x, y, z) or (value, x, y, z, photonID) if photon IDs were logged."""
        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, transform: TransformFn = _noTransform
    ) -> Optional[np.ndarray]:
        if key and key.solidLabel:
            if not self._keyExists(key):
                return None
            container: ListArrayContainer = getattr(self._data[key], dataType.value)
            return transform(key, container.getData())
        else:
            container = ListArrayContainer()
            for interactionData in self._data.values():
                points: Optional[ListArrayContainer] = getattr(interactionData, dataType.value)
                if points is None:
                    continue
                container.extend(points)
            if len(container) == 0:
                return None
            return container.getData()

    def _keyExists(self, key: InteractionKey) -> bool:
        if key.solidLabel not in self.getStoredSolidLabels():
            warnings.warn(
                f"No data stored for solid labeled '{key.solidLabel}'. Available: {self.getStoredSolidLabels()}. "
            )
        elif key.surfaceLabel and key.surfaceLabel not in self.getStoredSurfaceLabels(key.solidLabel):
            warnings.warn(
                f"No data stored for surface labeled '{key.surfaceLabel}' for solid '{key.solidLabel}'. "
                f"Available: {self.getStoredSurfaceLabels(key.solidLabel)}. "
            )
        if key in self._data:
            return True
        return False

    def save(self, filepath: str = None):
        if filepath is None and self._filepath is None:
            filepath = self.DEFAULT_LOGGER_PATH
            warnings.warn(f"No filepath specified. Saving to {filepath}.")
        elif filepath is None:
            filepath = self._filepath

        with open(filepath, "wb") as file:
            pickle.dump((self._data, self.info, self._labels), file)

    def load(self, filepath: str):
        self._filepath = filepath

        if not os.path.exists(filepath):
            warnings.warn(
                "No logger file found at '{}'. No data loaded, but it will create a new file "
                "at this location if the logger is saved later on.".format(filepath)
            )
            return

        with open(filepath, "rb") as file:
            self._data, self.info, self._labels = pickle.load(file)

    @property
    def hasFilePath(self):
        return self._filepath is not None

    @property
    def nDataPoints(self) -> int:
        return sum(len(data.dataPoints) for data in self._data.values())

getRawDataPoints(key=None)

All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) or (n, 5) where the second axis is (value, x, y, z) or (value, x, y, z, photonID) if photon IDs were logged.

Source code in pytissueoptics/scene/logger/logger.py

def getRawDataPoints(self, key: InteractionKey = None) -> np.ndarray:
    """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) or (n, 5) where
    the second axis is (value, x, y, z) or (value, x, y, z, photonID) if photon IDs were logged."""
    return self._getData(DataType.DATA_POINT, key)

getSeenSolidLabels()

Returns a list of all solid labels that have been logged in the past even if the data was discarded.

Source code in pytissueoptics/scene/logger/logger.py

def getSeenSolidLabels(self) -> List[str]:
    """Returns a list of all solid labels that have been logged in the past
    even if the data was discarded."""
    return list(self._labels.keys())

getSeenSurfaceLabels(solidLabel)

Returns a list of all surface labels that have been logged in the past for the given solid even if the data was discarded.

Source code in pytissueoptics/scene/logger/logger.py

def getSeenSurfaceLabels(self, solidLabel: str) -> List[str]:
    """Returns a list of all surface labels that have been logged in the past
    for the given solid even if the data was discarded."""
    return self._labels[solidLabel]

getStoredSolidLabels()

Returns a list of all solid labels that are currently stored in the logger.

Source code in pytissueoptics/scene/logger/logger.py

def getStoredSolidLabels(self) -> List[str]:
    """Returns a list of all solid labels that are currently stored in the logger."""
    return list(set(key.solidLabel for key in self._data.keys()))

getStoredSurfaceLabels(solidLabel)

Returns a list of all surface labels that are currently stored in the logger.

Source code in pytissueoptics/scene/logger/logger.py

def getStoredSurfaceLabels(self, solidLabel: str) -> List[str]:
    """Returns a list of all surface labels that are currently stored in the logger."""
    return [
        key.surfaceLabel
        for key in self._data.keys()
        if key.solidLabel == solidLabel and key.surfaceLabel is not None
    ]

logDataPointArray(array, key)

'array' must be of shape (n, 4) or (n, 5) where the second axis is (value, x, y, z) or (value, x, y, z, photonID). The photonID column is optional.

Source code in pytissueoptics/scene/logger/logger.py

def logDataPointArray(self, array: np.ndarray, key: InteractionKey):
    """'array' must be of shape (n, 4) or (n, 5) where the second axis is (value, x, y, z)
    or (value, x, y, z, photonID). The photonID column is optional."""
    assert array.shape[1] in [4, 5] and array.ndim == 2, "Data point array must be of shape (n, 4) or (n, 5)"
    self._appendData(array, DataType.DATA_POINT, key)

logPointArray(array, key=None)

'array' must be of shape (n, 3) where the second axis is (x, y, z)

Source code in pytissueoptics/scene/logger/logger.py

def logPointArray(self, array: np.ndarray, key: InteractionKey = None):
    """'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)

logSegmentArray(array, key=None)

'array' must be of shape (n, 6) where the second axis is (x1, y1, z1, x2, y2, z2)

Source code in pytissueoptics/scene/logger/logger.py

def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None):
    """'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)

Scene

Bases: Displayable

Source code in pytissueoptics/scene/scene/scene.py

class Scene(Displayable):
    def __init__(self, solids: List[Solid] = None, ignoreIntersections=False, worldMaterial=None):
        self._solids: List[Solid] = []
        self._ignoreIntersections = ignoreIntersections
        self._solidsContainedIn: Dict[str, List[str]] = {}
        self._worldMaterial = worldMaterial

        if solids:
            for solid in solids:
                self.add(solid)

        self.resetOutsideMaterial()

    def add(self, solid: Solid, position: Vector = None):
        if position:
            solid.translateTo(position)
        self._validateLabel(solid)
        if not self._ignoreIntersections:
            self._validatePosition(solid)
        self._solids.append(solid)

    @property
    def solids(self):
        return self._solids

    def addToViewer(self, viewer: Abstract3DViewer, representation="surface", colormap="bone", opacity=0.1, **kwargs):
        viewer.add(*self.solids, representation=representation, colormap=colormap, opacity=opacity, **kwargs)

    def getWorldEnvironment(self) -> Environment:
        return Environment(self._worldMaterial)

    def _validatePosition(self, newSolid: Solid):
        """Assert newSolid position is valid and make proper adjustments so that the
        material at each solid interface is well-defined."""
        if len(self._solids) == 0:
            return

        intersectingSuspects = self._findIntersectingSuspectsFor(newSolid)
        if len(intersectingSuspects) == 0:
            return

        intersectingSuspects.sort(key=lambda s: s.getBoundingBox().xMax - s.getBoundingBox().xMin, reverse=True)

        for otherSolid in intersectingSuspects:
            if newSolid.contains(*otherSolid.getVertices()):
                self._processContainedSolid(otherSolid, container=newSolid)
                break
            elif otherSolid.contains(*newSolid.getVertices()):
                self._processContainedSolid(newSolid, container=otherSolid)
            else:
                raise NotImplementedError(
                    "Cannot place a solid that partially intersects with an existing solid. "
                    "Since this might be underestimating containment, you can also create a "
                    "scene with 'ignoreIntersections=True' to ignore this error and manually "
                    "handle environments of contained solids with "
                    "containedSolid.setOutsideEnvironment(containerSolid.getEnvironment())."
                )

    def _processContainedSolid(self, solid: Solid, container: Solid):
        if container.isStack():
            containerEnv = self._getEnvironmentOfStackAt(solid.position, container)
        else:
            containerEnv = container.getEnvironment()
        solid.setOutsideEnvironment(containerEnv)

        containerLabel = containerEnv.solid.getLabel()
        if containerLabel not in self._solidsContainedIn:
            self._solidsContainedIn[containerLabel] = [solid.getLabel()]
        else:
            self._solidsContainedIn[containerLabel].append(solid.getLabel())

    def _validateLabel(self, solid):
        labelSet = set(s.getLabel() for s in self.solids)
        if solid.getLabel() not in labelSet:
            return
        idx = 2
        while f"{solid.getLabel()}_{idx}" in labelSet:
            idx += 1
        warnings.warn(
            f"A solid with label '{solid.getLabel()}' already exists in the scene. "
            f"Renaming to '{solid.getLabel()}_{idx}'."
        )
        solid.setLabel(f"{solid.getLabel()}_{idx}")

    def _findIntersectingSuspectsFor(self, solid) -> List[Solid]:
        solidBBox = solid.getBoundingBox()
        intersectingSuspects = []
        for otherSolid in self._solids:
            if solidBBox.intersects(otherSolid.getBoundingBox()):
                intersectingSuspects.append(otherSolid)
        return intersectingSuspects

    def getSolids(self) -> List[Solid]:
        return self._solids

    def getSolid(self, solidLabel: str) -> Solid:
        for solid in self._solids:
            if solid.getLabel().lower() == solidLabel.lower():
                return solid
            if not solid.isStack():
                continue
            for layerLabel in solid.getLayerLabels():
                if layerLabel.lower() == solidLabel.lower():
                    return solid
        raise ValueError(f"Solid '{solidLabel}' not found in scene. Available solids: {self.getSolidLabels()}")

    def getSolidLabels(self) -> List[str]:
        labels = []
        for solid in self._solids:
            if solid.isStack():
                labels.extend(solid.getLayerLabels())
            else:
                labels.append(solid.getLabel())
        return labels

    def getSurfaceLabels(self, solidLabel) -> List[str]:
        solid = self.getSolid(solidLabel)
        if solid is None:
            return []
        if solid.isStack() and solid.getLabel() != solidLabel:
            return solid.getLayerSurfaceLabels(solidLabel)
        return solid.surfaceLabels

    def getContainedSolidLabels(self, solidLabel: str) -> List[str]:
        return self._solidsContainedIn.get(solidLabel, [])

    def getPolygons(self) -> List[Polygon]:
        polygons = []
        for solid in self._solids:
            polygons.extend(solid.getPolygons())
        return polygons

    def getMaterials(self) -> list:
        materials = [self._worldMaterial]
        for solid in self._solids:
            if solid.isDetector:
                continue
            surfaceLabels = solid.surfaceLabels
            for surfaceLabel in surfaceLabels:
                material = solid.getPolygons(surfaceLabel)[0].insideEnvironment.material
                if material not in materials:
                    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

        bbox = self._solids[0].getBoundingBox().copy()
        for solid in self._solids[1:]:
            bbox.extendTo(solid.getBoundingBox())
        return bbox

    def resetOutsideMaterial(self):
        outsideEnvironment = self.getWorldEnvironment()
        for solid in self._solids:
            if self._isHidden(solid.getLabel()):
                continue
            solid.setOutsideEnvironment(outsideEnvironment)

    def _isHidden(self, solidLabel: str) -> bool:
        for hiddenLabels in self._solidsContainedIn.values():
            if solidLabel in hiddenLabels:
                return True
        return False

    def getEnvironmentAt(self, position: Vector) -> Environment:
        # First, recursively look if position is in a contained solid.
        for containerLabel in self._solidsContainedIn.keys():
            env = self._getEnvironmentOfContainerAt(position, containerLabel)
            if env is not None:
                return env

        for solid in self._solids:
            if solid.contains(position):
                if solid.isStack():
                    return self._getEnvironmentOfStackAt(position, solid)
                return solid.getEnvironment()
        return self.getWorldEnvironment()

    def _getEnvironmentOfContainerAt(self, position: Vector, containerLabel: str) -> Optional[Environment]:
        containerSolid = self.getSolid(containerLabel)
        if not containerSolid.contains(position):
            return None
        for containedLabel in self.getContainedSolidLabels(containerLabel):
            containedEnv = self._getEnvironmentOfContainerAt(position, containedLabel)
            if containedEnv:
                return containedEnv
        if containerSolid.isStack():
            return self._getEnvironmentOfStackAt(position, containerSolid)
        return containerSolid.getEnvironment()

    @staticmethod
    def _getEnvironmentOfStackAt(position: Vector, stack: Solid) -> Environment:
        """Returns the environment of the stack at the given position.

        To do that we first find the interface in the stack that is closest to the given position.
        At the same time we find on which side of the interface we are and return the environment
        of this side from any surface polygon.
        """
        environment = None
        closestDistance = sys.maxsize
        for surfaceLabel in stack.surfaceLabels:
            if INTERFACE_KEY not in surfaceLabel:
                continue
            planePolygon = stack.surfaces.getPolygons(surfaceLabel)[0]
            planeNormal = planePolygon.normal
            planePoint = planePolygon.vertices[0]
            v = position - planePoint
            distanceFromPlane = v.dot(planeNormal)
            if abs(distanceFromPlane) < closestDistance:
                closestDistance = abs(distanceFromPlane)
                isInside = distanceFromPlane < 0
                environment = planePolygon.insideEnvironment if isInside else planePolygon.outsideEnvironment
        return environment

    def __hash__(self):
        solidHash = hash(tuple(sorted([hash(s) for s in self._solids])))
        worldMaterialHash = hash(self._worldMaterial) if self._worldMaterial else 0
        return hash((solidHash, worldMaterialHash))

Sphere

Bases: Ellipsoid

The Sphere is the 3D analog to the circle. Meshing a sphere requires an infinite number of vertices. The position refers to the vector from global origin to its centroid. The radius of the sphere will determine the outermost distance from its centroid.

This class offers two possible methods to generate the sphere mesh. - With Quads: Specify the number of separation lines on the vertical axis and the horizontal axis of the sphere. - With Triangle: Specify the order of splitting. This will generate what is known as an IcoSphere.

Source code in pytissueoptics/scene/solids/sphere.py

class Sphere(Ellipsoid):
    """
    The Sphere is the 3D analog to the circle. Meshing a sphere requires an infinite number of vertices.
    The position refers to the vector from global origin to its centroid.
    The radius of the sphere will determine the outermost distance from its centroid.

    This class offers two possible methods to generate the sphere mesh.
    - With Quads: Specify the number of separation lines on the vertical axis and the horizontal axis of the sphere.
    - With Triangle: Specify the order of splitting. This will generate what is known as an IcoSphere.
    """

    def __init__(
        self,
        radius: float = 1.0,
        order: int = 3,
        position: Vector = Vector(0, 0, 0),
        material=None,
        label: str = "sphere",
        primitive: str = primitives.DEFAULT,
        smooth: bool = True,
    ):
        self._radius = radius

        super().__init__(
            a=radius,
            b=radius,
            c=radius,
            order=order,
            position=position,
            material=material,
            label=label,
            primitive=primitive,
            smooth=smooth,
        )

    @property
    def radius(self):
        return self._radius

    def _computeQuadMesh(self):
        raise NotImplementedError

    def contains(self, *vertices: Vector) -> bool:
        """Only returns true if all vertices are inside the minimum radius of the sphere
        (more restrictive with low order spheres)."""
        minRadius = self._getMinimumRadius()
        for vertex in vertices:
            relativeVertex = vertex - self.position
            if relativeVertex.getNorm() >= minRadius:
                return False
        return True

    def _getMinimumRadius(self) -> float:
        return (1 - self._getRadiusError()) * self._radius

    def _radiusTowards(self, vertex) -> float:
        return self.radius

    def _geometryParams(self) -> dict:
        return {
            "radius": self._radius,
            "order": self._order,
        }

contains(*vertices)

Only returns true if all vertices are inside the minimum radius of the sphere (more restrictive with low order spheres).

Source code in pytissueoptics/scene/solids/sphere.py

def contains(self, *vertices: Vector) -> bool:
    """Only returns true if all vertices are inside the minimum radius of the sphere
    (more restrictive with low order spheres)."""
    minRadius = self._getMinimumRadius()
    for vertex in vertices:
        relativeVertex = vertex - self.position
        if relativeVertex.getNorm() >= minRadius:
            return False
    return True

SymmetricLens

Bases: ThickLens

A symmetrical thick lens of focal length f in air.

Source code in pytissueoptics/scene/solids/lens.py

class SymmetricLens(ThickLens):
    """A symmetrical thick lens of focal length `f` in air."""

    def __init__(
        self,
        f: float,
        diameter: float,
        thickness: float,
        material: RefractiveMaterial,
        position: Vector = Vector(0, 0, 0),
        label: str = "lens",
        primitive: str = primitives.DEFAULT,
        smooth: bool = True,
        u: int = 24,
        v: int = 2,
        s: int = 24,
    ):
        # For thick lenses, the focal length is given by the lensmaker's equation:
        # 1/f = (n - 1) * (1/R1 - 1/R2 + (n - 1) * d / (n * R1 * R2))
        # with R2 = -R1, we get the following quadratic equation to solve:
        # 1/f = (n - 1) * (2/R - (n - 1) * d / (n * R^2))
        n = material.n
        p = math.sqrt(f * n * (f * n - thickness)) * (n - 1) / n
        R = f * (n - 1) + p * np.sign(f)
        super().__init__(R, -R, diameter, thickness, position, material, label, primitive, smooth, u, v, s)

ThickLens

Bases: Cylinder

The Lens is defined by a front radius, a back radius, a diameter and a thickness (along its center). A positive frontRadius means that the front surface is convex. This is reversed for the back surface. A flat surface is obtained by setting the corresponding radius to 0 or math.inf. The position refers to the vector from global origin to its centroid. The generated mesh will be divided into the following subgroups: Front, Side and Back. By default, front will point towards the negative z-axis.

Source code in pytissueoptics/scene/solids/lens.py

class ThickLens(Cylinder):
    """
    The Lens is defined by a front radius, a back radius, a diameter and a thickness (along its center).
    A positive frontRadius means that the front surface is convex. This is reversed for the back surface.
    A flat surface is obtained by setting the corresponding radius to 0 or math.inf.
    The position refers to the vector from global origin to its centroid.
    The generated mesh will be divided into the following subgroups: Front, Side and Back.
    By default, front will point towards the negative z-axis.
    """

    def __init__(
        self,
        frontRadius: float,
        backRadius: float,
        diameter: float,
        thickness: float,
        position: Vector = Vector(0, 0, 0),
        material=None,
        label: str = "thick lens",
        primitive: str = primitives.DEFAULT,
        smooth: bool = True,
        u: int = 24,
        v: int = 2,
        s: int = 24,
    ):
        frontRadius = frontRadius if frontRadius != 0 else math.inf
        backRadius = backRadius if backRadius != 0 else math.inf
        if abs(frontRadius) <= diameter / 2:
            raise ValueError(
                f"Front radius must be greater than the lens radius. Front radius: {frontRadius}, "
                f"lens radius: {diameter / 2}"
            )
        if abs(backRadius) <= diameter / 2:
            raise ValueError(
                f"Back radius must be greater than the lens radius. Back radius: {backRadius}, "
                f"lens radius: {diameter / 2}"
            )

        self._diameter = diameter
        self._frontRadius = frontRadius
        self._backRadius = backRadius
        self._thickness = thickness

        length = self._computeEdgeThickness(thickness)
        super().__init__(
            radius=diameter / 2,
            length=length,
            u=u,
            v=v,
            s=s,
            position=position,
            material=material,
            label=label,
            primitive=primitive,
            smooth=smooth,
        )

    @property
    def _hasFrontCurvature(self) -> bool:
        return self._frontRadius != 0 and self._frontRadius != math.inf

    @property
    def _hasBackCurvature(self) -> bool:
        return self._backRadius != 0 and self._backRadius != math.inf

    def _computeEdgeThickness(self, centerThickness) -> float:
        """Returns the thickness of the lens on its side. This is the length required to build the base cylinder
        before applying surface curvature."""
        dt1, dt2 = 0, 0
        if self._hasFrontCurvature:
            dt1 = abs(self._frontRadius) - math.sqrt(self._frontRadius**2 - self._diameter**2 / 4)
            dt1 *= np.sign(self._frontRadius)
        if self._hasBackCurvature:
            dt2 = abs(self._backRadius) - math.sqrt(self._backRadius**2 - self._diameter**2 / 4)
            dt2 *= -np.sign(self._backRadius)
        edgeThickness = centerThickness - dt1 - dt2
        if edgeThickness < 0:
            raise ValueError("Desired center thickness is too small for the given radii and diameter.")
        return edgeThickness

    @property
    def centerThickness(self) -> float:
        return (self._frontCenter - self._backCenter).getNorm()

    @property
    def edgeThickness(self) -> float:
        return self._length

    @property
    def frontRadius(self) -> float:
        return self._frontRadius

    @property
    def backRadius(self) -> float:
        return self._backRadius

    @property
    def focalLength(self) -> float:
        """Returns the focal length of the lens in air. Requires a refractive material to be defined."""
        if self._material is None or not issubclass(type(self._material), RefractiveMaterial):
            raise ValueError("Cannot compute focal length without refractive material defined.")
        # For thick lenses, the focal length is given by the lensmaker's equation:
        # 1/f = (n - 1) * (1/R1 - 1/R2 + (n - 1) * d / (n * R1 * R2))
        n = self._material.n
        R1 = self._frontRadius
        R2 = self._backRadius
        d = self.centerThickness
        if n == 1:
            return math.inf
        if self._hasFrontCurvature and self._hasBackCurvature:
            return 1 / ((n - 1) * (1 / R1 - 1 / R2 + (n - 1) * d / (n * R1 * R2)))
        if self._hasFrontCurvature:
            return R1 / (n - 1)
        if self._hasBackCurvature:
            return R2 / (n - 1)
        return math.inf

    def _computeVerticesOfLayers(self) -> tuple:
        frontLayers, lateralLayers, backLayers = super()._computeVerticesOfLayers()

        if self._hasFrontCurvature:
            frontVertices = list(itertools.chain.from_iterable(frontLayers)) + [self._frontCenter]
            self._applyCurvature(self._frontRadius, frontVertices)
        if self._hasBackCurvature:
            backVertices = list(itertools.chain.from_iterable(backLayers)) + [self._backCenter]
            self._applyCurvature(self._backRadius, backVertices)

        if self._frontCenter.z > self._backCenter.z:
            raise ValueError("Not a valid lens: curved surfaces intersect.")
        return frontLayers, lateralLayers, backLayers

    def _applyCurvature(self, radius: float, vertices: List[Vertex]):
        # At this point, all vertices are on the same z plane.
        surfaceZ = vertices[0].z
        # The sphere origin is simply found by setting the z coordinate so that the distance to
        # the surface perimeter equals the desired radius.
        sphereOrigin = Vector(0, 0, surfaceZ + math.sqrt(radius**2 - self._radius**2) * np.sign(radius))
        for vertex in vertices:
            direction = vertex - sphereOrigin
            direction.normalize()
            vertex.update(*(sphereOrigin + direction * abs(radius)).array)

    def smooth(self, surfaceLabel: str = None, reset: bool = True):
        if surfaceLabel:
            return super(Cylinder, self).smooth(surfaceLabel, reset)
        for surfaceLabel in ["front", "back"]:
            self.smooth(surfaceLabel, reset=False)

    def _geometryParams(self) -> dict:
        return {
            "diameter": self._diameter,
            "frontRadius": self._frontRadius,
            "backRadius": self._backRadius,
            "thickness": self._thickness,
            "length": self._length,
            "u": self._u,
            "v": self._v,
            "s": self._s,
        }

focalLength property

Returns the focal length of the lens in air. Requires a refractive material to be defined.

Vector

Basic implementation of a mutable 3D Vector. It implements most of the basic vector operation. Mutability is necessary when working with shared object references for expected behavior.

Source code in pytissueoptics/scene/geometry/vector.py

class Vector:
    """
    Basic implementation of a mutable 3D Vector. It implements most of the basic vector operation.
    Mutability is necessary when working with shared object references for expected behavior.
    """

    def __init__(self, x: float = 0, y: float = 0, z: float = 0):
        self._x = x
        self._y = y
        self._z = z

    @property
    def x(self):
        return self._x

    @property
    def y(self):
        return self._y

    @property
    def z(self):
        return self._z

    def __repr__(self):
        return f"<Vector>:({self._x}, {self._y}, {self._z})"

    def __iter__(self):
        return iter((self._x, self._y, self._z))

    def __eq__(self, other: "Vector"):
        tol = 1e-5
        if (
            math.isclose(other._x, self._x, abs_tol=tol)
            and math.isclose(other._y, self._y, abs_tol=tol)
            and math.isclose(other._z, self._z, abs_tol=tol)
        ):
            return True
        else:
            return False

    def __sub__(self, other: "Vector") -> "Vector":
        return Vector(self._x - other._x, self._y - other._y, self._z - other._z)

    def __add__(self, other: "Vector") -> "Vector":
        return Vector(self._x + other._x, self._y + other._y, self._z + other._z)

    def __mul__(self, scalar: float) -> "Vector":
        return Vector(self._x * scalar, self._y * scalar, self._z * scalar)

    def __truediv__(self, scalar: float) -> "Vector":
        return Vector(self._x / scalar, self._y / scalar, self._z / scalar)

    def add(self, other: "Vector"):
        self._x += other._x
        self._y += other._y
        self._z += other._z

    def subtract(self, other: "Vector"):
        self._x -= other._x
        self._y -= other._y
        self._z -= other._z

    def multiply(self, scalar: float):
        self._x *= scalar
        self._y *= scalar
        self._z *= scalar

    def divide(self, scalar: float):
        self._x /= scalar
        self._y /= scalar
        self._z /= scalar

    def getNorm(self) -> float:
        return (self._x**2 + self._y**2 + self._z**2) ** (1 / 2)

    def normalize(self):
        norm = self.getNorm()
        if norm != 0:
            self._x = self._x / norm
            self._y = self._y / norm
            self._z = self._z / norm

    def cross(self, other: "Vector") -> "Vector":
        ux, uy, uz = self._x, self._y, self._z
        vx, vy, vz = other._x, other._y, other._z
        return Vector(uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx)

    def dot(self, other: "Vector") -> float:
        return self._x * other._x + self._y * other._y + self._z * other._z

    @property
    def array(self) -> list:
        return [self._x, self._y, self._z]

    def update(self, x: float, y: float, z: float):
        self._x = x
        self._y = y
        self._z = z

    def copy(self) -> "Vector":
        return Vector(self._x, self._y, self._z)

    def rotateAround(self, unitAxis: "Vector", theta: float):
        """
        Rotate the vector around `unitAxis` by `theta` radians. Assumes the axis to be a unit vector.
        Uses Rodrigues' rotation formula.
        """
        # This is the most expensive (and most common)
        # operation when performing Monte Carlo in tissue
        # (15% of time spent here). It is difficult to optimize without
        # making it even less readable than it currently is
        # http://en.wikipedia.org/wiki/Rotation_matrix
        #
        # Several options were tried in the past such as
        # external not-so-portable C library, unreadable
        # shortcuts, sine and cosine lookup tables, etc...
        # and the performance gain was minimal (<20%).
        # For now, this is the best, most readable solution.

        cost = math.cos(theta)
        sint = math.sin(theta)
        one_cost = 1 - cost

        ux = unitAxis.x
        uy = unitAxis.y
        uz = unitAxis.z

        X = self._x
        Y = self._y
        Z = self._z

        x = (
            (cost + ux * ux * one_cost) * X
            + (ux * uy * one_cost - uz * sint) * Y
            + (ux * uz * one_cost + uy * sint) * Z
        )
        y = (
            (uy * ux * one_cost + uz * sint) * X
            + (cost + uy * uy * one_cost) * Y
            + (uy * uz * one_cost - ux * sint) * Z
        )
        z = (
            (uz * ux * one_cost - uy * sint) * X
            + (uz * uy * one_cost + ux * sint) * Y
            + (cost + uz * uz * one_cost) * Z
        )

        self.update(x, y, z)

    def getAnyOrthogonal(self) -> "Vector":
        if abs(self._z) < abs(self._x):
            return Vector(self._y, -self._x, 0)

        return Vector(0, -self._z, self._y)

    def __hash__(self):
        return hash((self._x, self._y, self._z))

rotateAround(unitAxis, theta)

Rotate the vector around unitAxis by theta radians. Assumes the axis to be a unit vector. Uses Rodrigues' rotation formula.

Source code in pytissueoptics/scene/geometry/vector.py

def rotateAround(self, unitAxis: "Vector", theta: float):
    """
    Rotate the vector around `unitAxis` by `theta` radians. Assumes the axis to be a unit vector.
    Uses Rodrigues' rotation formula.
    """
    # This is the most expensive (and most common)
    # operation when performing Monte Carlo in tissue
    # (15% of time spent here). It is difficult to optimize without
    # making it even less readable than it currently is
    # http://en.wikipedia.org/wiki/Rotation_matrix
    #
    # Several options were tried in the past such as
    # external not-so-portable C library, unreadable
    # shortcuts, sine and cosine lookup tables, etc...
    # and the performance gain was minimal (<20%).
    # For now, this is the best, most readable solution.

    cost = math.cos(theta)
    sint = math.sin(theta)
    one_cost = 1 - cost

    ux = unitAxis.x
    uy = unitAxis.y
    uz = unitAxis.z

    X = self._x
    Y = self._y
    Z = self._z

    x = (
        (cost + ux * ux * one_cost) * X
        + (ux * uy * one_cost - uz * sint) * Y
        + (ux * uz * one_cost + uy * sint) * Z
    )
    y = (
        (uy * ux * one_cost + uz * sint) * X
        + (cost + uy * uy * one_cost) * Y
        + (uy * uz * one_cost - ux * sint) * Z
    )
    z = (
        (uz * ux * one_cost - uy * sint) * X
        + (uz * uy * one_cost + ux * sint) * Y
        + (cost + uz * uz * one_cost) * Z
    )

    self.update(x, y, z)