Source code for eradiate.scenes.biosphere._leaf_cloud

from __future__ import annotations

import os
import warnings

import attrs
import mitsuba as mi
import numpy as np
import pint
import pinttr
import scipy as sp
import scipy.special

from ._core import CanopyElement
from ..core import SceneElement, traverse
from ..spectra import Spectrum, spectrum_factory
from ... import validators
from ...attrs import documented, get_doc, parse_docs
from ...kernel import TypeIdLookupStrategy, UpdateParameter
from ...units import unit_context_config as ucc
from ...units import unit_context_kernel as uck
from ...units import unit_registry as ureg


def _sample_lad(mu, nu, rng):
    """
    Generate an angle sample from the Leaf angle distribution function according to
    :cite:`GoelStrebel1984`, using the rejection method.
    """

    while True:
        rands = rng.random(2)
        theta_candidate = rands[0] * np.pi / 2.0
        gs_lad = (
            2.0
            / np.pi
            * sp.special.gamma(mu + nu)
            / (sp.special.gamma(mu) * sp.special.gamma(mu))
            * pow((1 - (2 * theta_candidate) / np.pi), mu - 1)
            * pow((2 * theta_candidate) / np.pi, nu - 1)
        )

        # scaling factor for the rejection method set to 2.0 to encompass the
        # entire distribution
        if rands[1] * 2.0 <= gs_lad:
            return theta_candidate


@ureg.wraps(ureg.m, (None, ureg.m, ureg.m, None))
def _leaf_cloud_positions_cuboid(n_leaves, l_horizontal, l_vertical, rng):
    """
    Compute leaf positions for a cuboid-shaped leaf cloud (square footprint).
    """
    positions = np.empty((n_leaves, 3))

    for i in range(n_leaves):
        rand = rng.random(3)
        positions[i, :] = [
            rand[0] * l_horizontal - 0.5 * l_horizontal,
            rand[1] * l_horizontal - 0.5 * l_horizontal,
            rand[2] * l_vertical,
        ]

    return positions


@ureg.wraps(ureg.m, (None, ureg.m, ureg.m, ureg.m, None, None))
def _leaf_cloud_positions_cuboid_avoid_overlap(
    n_leaves, l_horizontal, l_vertical, leaf_radius, n_attempts, rng
):
    """
    Compute leaf positions for a cuboid-shaped leaf cloud (square footprint).
    This function also performs conservative collision checks to avoid leaf
    overlapping. This process might take a very long time, if the parameters
    specify a very dense leaf cloud. Consider using
    :func:`_leaf_cloud_positions_cuboid`.
    """
    try:
        import aabbtree
    except ModuleNotFoundError:
        warnings.warn(
            "To use the collision detection feature, you must install AABBTree.\n"
            "See instructions on https://aabbtree.readthedocs.io/#installation."
        )
        raise

    n_attempts = int(n_attempts)  # For safety, ensure conversion to int

    # try placing the leaves such that they do not overlap by creating
    # axis-aligned bounding boxes and checking them for intersection
    positions = np.empty((n_leaves, 3))
    tree = aabbtree.AABBTree()

    for i in range(n_leaves):
        for j in range(n_attempts):
            rand = rng.random(3)
            pos_candidate = [
                rand[0] * l_horizontal - 0.5 * l_horizontal,
                rand[1] * l_horizontal - 0.5 * l_horizontal,
                rand[2] * l_vertical,
            ]
            aabb = aabbtree.AABB(
                [
                    (pos_candidate[0] - leaf_radius, pos_candidate[0] + leaf_radius),
                    (pos_candidate[1] - leaf_radius, pos_candidate[1] + leaf_radius),
                    (pos_candidate[2] - leaf_radius, pos_candidate[2] + leaf_radius),
                ]
            )
            if i == 0:
                positions[i, :] = pos_candidate
                tree.add(aabb)
                break
            else:
                if not tree.does_overlap(aabb):
                    positions[i, :] = pos_candidate
                    tree.add(aabb)
                    break
        else:
            raise RuntimeError(
                "unable to place all leaves: the specified canopy might be too dense"
            )

    return positions


@ureg.wraps(ureg.m, (None, None, ureg.m, ureg.m, ureg.m))
def _leaf_cloud_positions_ellipsoid(n_leaves: int, rng, a: float, b: float, c: float):
    """
    Compute leaf positions for an ellipsoid leaf cloud.
    The ellipsoid follows the equation:
    :math:`\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2}= 1`
    """

    positions = []

    while len(positions) < n_leaves:
        rand = rng.random(3)
        x = (rand[0] - 0.5) * 2 * a
        y = (rand[1] - 0.5) * 2 * b
        z = (rand[2] - 0.5) * 2 * c

        if (x**2 / a**2) + (y**2 / b**2) + (z**2 / c**2) <= 1.0:
            positions.append([x, y, z])

    return positions


@ureg.wraps(ureg.m, (None, ureg.m, ureg.m, None))
def _leaf_cloud_positions_cylinder(n_leaves, radius, l_vertical, rng):
    """
    Compute leaf positions for a cylinder-shaped leaf cloud (vertical
    orientation).
    """

    positions = np.empty((n_leaves, 3))

    for i in range(n_leaves):
        rand = rng.random(3)
        phi = rand[0] * 2 * np.pi
        r = rand[1] * radius
        z = rand[2] * l_vertical
        positions[i, :] = [r * np.cos(phi), r * np.sin(phi), z]

    return positions


@ureg.wraps(ureg.m, (None, ureg.m, ureg.m, None))
def _leaf_cloud_positions_cone(n_leaves, radius, l_vertical, rng):
    """
    Compute leaf positions for a cone-shaped leaf cloud (vertical
    orientation, tip pointing towards positive z).
    """

    positions = np.empty((n_leaves, 3))

    # uniform cone sampling from here:
    # https://stackoverflow.com/questions/41749411/uniform-sampling-by-volume-within-a-cone
    for i in range(n_leaves):
        rand = rng.random(3)
        h = l_vertical * (rand[0] ** (1 / 3))
        r = radius / l_vertical * h * np.sqrt(rand[1])
        phi = rand[2] * 2 * np.pi
        positions[i, :] = [r * np.cos(phi), r * np.sin(phi), l_vertical - h]

    return positions


@ureg.wraps(None, (None, None, None, None))
def _leaf_cloud_orientations(n_leaves, mu, nu, rng):
    """Compute leaf orientations."""
    orientations = np.empty((n_leaves, 3))
    for i in range(np.shape(orientations)[0]):
        theta = _sample_lad(mu, nu, rng)
        phi = rng.random() * 2.0 * np.pi

        orientations[i, :] = [
            np.sin(theta) * np.cos(phi),
            np.sin(theta) * np.sin(phi),
            np.cos(theta),
        ]

    return orientations


@ureg.wraps(ureg.m, (None, ureg.m))
def _leaf_cloud_radii(n_leaves, leaf_radius):
    """Compute leaf radii."""
    return np.full((n_leaves,), leaf_radius)


@parse_docs
@attrs.define
class LeafCloudParams:
    """
    Base class to implement advanced parameter checking for :class:`.LeafCloud`
    generators.
    """

    _id = documented(
        attrs.field(default="leaf_cloud"),
        doc="Leaf cloud identifier.",
        type="str",
        default='"leaf_cloud"',
    )

    _leaf_reflectance = documented(
        attrs.field(default=0.5), doc="Leaf reflectance.", type="float", default="0.5"
    )

    _leaf_transmittance = documented(
        attrs.field(default=0.5), doc="Leaf transmittance.", type="float", default="0.5"
    )

    _mu = documented(
        attrs.field(default=1.066),
        doc="First parameter of the inverse beta distribution approximation used "
        "to generate leaf orientations.",
        type="float",
        default="1.066",
    )

    _nu = documented(
        attrs.field(default=1.853),
        doc="Second parameter of the inverse beta distribution approximation used "
        "to generate leaf orientations.",
        type="float",
        default="1.853",
    )

    _n_leaves = documented(
        attrs.field(default=None), doc="Number of leaves.", type="int"
    )

    _leaf_radius = documented(
        pinttr.field(default=None, units=ucc.deferred("length")),
        doc="Leaf radius.\n\nUnit-enabled field (default: ucc['length']).",
        type="float",
    )

    def update(self):
        try:
            for field in [x.name.lstrip("_") for x in self.__attrs_attrs__]:
                self.__getattribute__(field)
        except Exception as e:
            raise Exception(
                f"cannot compute field '{field}', parameter set is likely under-constrained"
            ) from e

    def __attrs_post_init__(self):
        self.update()

    @property
    def id(self):
        return self._id

    @property
    def leaf_reflectance(self):
        return self._leaf_reflectance

    @property
    def leaf_transmittance(self):
        return self._leaf_transmittance

    @property
    def nu(self):
        return self._nu

    @property
    def mu(self):
        return self._mu

    @property
    def n_leaves(self):
        return self._n_leaves

    @property
    def leaf_radius(self):
        return self._leaf_radius


[docs]@parse_docs @attrs.define class CuboidLeafCloudParams(LeafCloudParams): """ Advanced parameter checking class for the cuboid :class:`.LeafCloud` generator. Some of the parameters can be inferred from each other. Parameters defined below can be used (without leading underscore) as keyword arguments to the :meth:`.LeafCloud.cuboid` class method constructor. Parameters without defaults are connected by a dependency graph used to compute required parameters (outlined in the figure below). The following parameter sets are valid: * ``n_leaves``, ``leaf_radius``, ``l_horizontal``, ``l_vertical``; * ``lai``, ``leaf_radius``, ``l_horizontal``, ``l_vertical``; * ``lai``, ``leaf_radius``, ``l_horizontal``, ``hdo``, ``hvr``; * and more! .. only:: latex .. figure:: ../../../../fig/cuboid_leaf_cloud_params.png .. only:: not latex .. figure:: ../../../../fig/cuboid_leaf_cloud_params.svg Warnings -------- In case of over-specification, no consistency check is performed. See Also -------- :meth:`.LeafCloud.cuboid` """ _l_horizontal = documented( pinttr.field(default=None, units=ucc.deferred("length")), doc="Leaf cloud horizontal extent. *Suggested default: 30 m.*\n" "\n" "Unit-enabled field (default: ucc['length']).", type="float", ) _l_vertical = documented( pinttr.field(default=None, units=ucc.deferred("length")), doc="Leaf cloud vertical extent. *Suggested default: 3 m.*\n" "\n" "Unit-enabled field (default: ucc['length']).", type="float", ) _lai = documented( pinttr.field(default=None, units=ureg.dimensionless), doc="Leaf cloud leaf area index (LAI). *Physical range: [0, 10]; " "suggested default: 3.*\n" "\n" "Unit-enabled field (default: ucc['dimensionless']).", type="float", ) _hdo = documented( pinttr.field(default=None, units=ucc.deferred("length")), doc="Mean horizontal distance between leaves.\n" "\n" "Unit-enabled field (default: ucc['length']).", type="float", ) _hvr = documented( pinttr.field(default=None), doc="Ratio of mean horizontal leaf distance and vertical leaf cloud extent. " "*Suggested default: 0.1.*", type="float", ) @property def n_leaves(self): if self._n_leaves is None: self._n_leaves = int( self.lai * (self.l_horizontal / self.leaf_radius) ** 2 / np.pi ) return self._n_leaves @property def lai(self): if self._lai is None: self._lai = ( np.pi * (self.leaf_radius / self.l_horizontal) ** 2 * self.n_leaves ) return self._lai @property def leaf_radius(self): if self._leaf_radius is None: self._leaf_radius = ( np.sqrt(self.lai / (self.n_leaves * np.pi)) * self.l_horizontal ) return self._leaf_radius @property def l_horizontal(self): if self._l_horizontal is None: self._l_horizontal = ( np.pi * self.leaf_radius**2 * self.n_leaves / self.lai ) return self._l_horizontal @property def l_vertical(self): if self._l_vertical is None: self._l_vertical = ( self.lai * self.hdo**3 / (np.pi * self.leaf_radius**2 * self.hvr) ) return self._l_vertical @property def hdo(self): return self._hdo @property def hvr(self): return self._hvr def __str__(self): result = [] for field in [ "id", "lai", "leaf_radius", "l_horizontal", "l_vertical", "n_leaves", "leaf_reflectance", "leaf_transmittance", ]: value = self.__getattribute__(field) result.append(f"{field}={value.__repr__()}") return f"CuboidLeafCloudParams({', '.join(result)})"
[docs]@parse_docs @attrs.define class SphereLeafCloudParams(LeafCloudParams): """ Advanced parameter checking class for the sphere :class:`.LeafCloud` generator. See Also -------- :meth:`.LeafCloud.sphere` """ _radius = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) @property def radius(self): return self._radius
[docs]@parse_docs @attrs.define class EllipsoidLeafCloudParams(LeafCloudParams): """ Advanced parameter checking class for the ellipsoid :class:`.LeafCloud` generator. Parameters ``a``, ``b`` and ``c`` denote the ellipsoid's half axes along the x, y, and z directions respectively. If either ``b`` or ``c`` are not set by the user, they default to being equal to ``a``. Accordingly, a sphere of radius ``r`` can be parametrized by setting ``a=r``. See Also -------- :meth:`.LeafCloud.ellipsoid` """ _a = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) _b = documented( pinttr.field(default=None, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) _c = documented( pinttr.field(default=None, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) @property def a(self): if self._a <= 0: raise ValueError( "Ellipsoid half axis parameters must be strictly larger than zero!" ) return self._a @property def b(self): if self._b is None: self._b = self.a elif self._b <= 0: raise ValueError( "Ellipsoid half axis parameters must be strictly larger than zero!" ) return self._b @property def c(self): if self._c is None: self._c = self.a elif self._c <= 0: raise ValueError( "Ellipsoid half axis parameters must be strictly larger than zero!" ) return self._c
[docs]@parse_docs @attrs.define class CylinderLeafCloudParams(LeafCloudParams): """ Advanced parameter checking class for the cylinder :class:`.LeafCloud` generator. See Also -------- :meth:`.LeafCloud.cylinder` """ _radius = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) _l_vertical = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud vertical extent.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) @property def radius(self): return self._radius @property def l_vertical(self): return self._l_vertical
[docs]@parse_docs @attrs.define class ConeLeafCloudParams(LeafCloudParams): """ Advanced parameter checking class for the cone :class:`.LeafCloud` generator. See Also -------- :meth:`.LeafCloud.cone` """ _radius = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud radius.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) _l_vertical = documented( pinttr.field(default=1.0 * ureg.m, units=ucc.deferred("length")), doc="Leaf cloud vertical extent.\n\nUnit-enabled field (default: ucc[length]).", type="float", default="1 m", ) @property def radius(self): return self._radius @property def l_vertical(self): return self._l_vertical
[docs]@parse_docs @attrs.define(eq=False, slots=False) class LeafCloud(CanopyElement): """ A container class for leaf clouds in abstract discrete canopies. Holds parameters completely characterizing the leaf cloud's leaves. In practice, this class should rarely be instantiated directly using its constructor. Instead, several class method constructors are available: * generators create leaf clouds from a set of parameters: * :meth:`.LeafCloud.cone`; * :meth:`.LeafCloud.cuboid`; * :meth:`.LeafCloud.cylinder`; * :meth:`.LeafCloud.ellipsoid`; * :meth:`.LeafCloud.sphere`; * :meth:`.LeafCloud.from_file` loads leaf positions and orientations from a text file. .. admonition:: Class method constructors .. autosummary:: cuboid cylinder ellipsoid from_file sphere """ # -------------------------------------------------------------------------- # Fields # -------------------------------------------------------------------------- id: str | None = documented( attrs.field( default="leaf_cloud", validator=attrs.validators.optional(attrs.validators.instance_of(str)), ), doc=get_doc(SceneElement, "id", "doc"), type=get_doc(SceneElement, "id", "type"), init_type=get_doc(SceneElement, "id", "init_type"), default="'leaf_cloud'", ) leaf_positions: pint.Quantity = documented( pinttr.field(factory=list, units=ucc.deferred("length")), doc="Leaf positions in cartesian coordinates as a (n, 3)-array.\n" "\n" "Unit-enabled field (default: ucc['length']).", type="quantity", init_type="array-like", default="[]", ) leaf_orientations: np.ndarray = documented( attrs.field(factory=list, converter=np.array), doc="Leaf orientations (normal vectors) in Cartesian coordinates as a " "(n, 3)-array.", type="ndarray", default="[]", ) leaf_radii: pint.Quantity = documented( pinttr.field( factory=list, validator=[ pinttr.validators.has_compatible_units, attrs.validators.deep_iterable(member_validator=validators.is_positive), ], units=ucc.deferred("length"), ), doc="Leaf radii as a n-array.\n\nUnit-enabled field (default: ucc[length]).", init_type="array-like", type="quantity", default="[]", ) @leaf_positions.validator @leaf_orientations.validator def _positions_orientations_validator(self, attribute, value): if not len(value): return if not value.ndim == 2 or value.shape[1] != 3: raise ValueError( f"While validating {attribute.name}: shape should be (N, 3), " f"got {value.shape}" ) @leaf_positions.validator @leaf_orientations.validator @leaf_radii.validator def _positions_orientations_radii_validator(self, attribute, value): if not ( len(self.leaf_positions) == len(self.leaf_orientations) == len(self.leaf_radii) ): raise ValueError( f"While validating {attribute.name}: " f"leaf_positions, leaf_orientations and leaf_radii must have the " f"same length. Got " f"len(leaf_positions) = {len(self.leaf_positions)}, " f"len(leaf_orientations) = {len(self.leaf_orientations)}, " f"len(leaf_radii) = {len(self.leaf_radii)}." ) leaf_reflectance: Spectrum = documented( attrs.field( default=0.5, converter=spectrum_factory.converter("reflectance"), validator=[ attrs.validators.instance_of(Spectrum), validators.has_quantity("reflectance"), ], ), doc="Reflectance spectrum of the leaves in the cloud. " "Must be a reflectance spectrum (dimensionless).", type=":class:`.Spectrum`", init_type=":class:`.Spectrum` or dict", default="0.5", ) leaf_transmittance: Spectrum = documented( attrs.field( default=0.5, converter=spectrum_factory.converter("transmittance"), validator=[ attrs.validators.instance_of(Spectrum), validators.has_quantity("transmittance"), ], ), doc="Transmittance spectrum of the leaves in the cloud. " "Must be a transmittance spectrum (dimensionless).", type=":class:`.Spectrum`", init_type=":class:`.Spectrum` or dict", default="0.5", ) # -------------------------------------------------------------------------- # Properties and accessors # --------------------------------------------------------------------------
[docs] def n_leaves(self) -> int: """ Returns ------- int Number of leaves in the leaf cloud. """ return len(self.leaf_positions)
[docs] def surface_area(self) -> pint.Quantity: """ Returns ------- quantity Total surface area. """ return np.sum(np.pi * self.leaf_radii * self.leaf_radii).squeeze()
# -------------------------------------------------------------------------- # Constructors # --------------------------------------------------------------------------
[docs] @classmethod def cuboid( cls, seed: int = 12345, avoid_overlap: bool = False, **kwargs ) -> LeafCloud: """ Generate a leaf cloud with an axis-aligned cuboid shape (and a square footprint on the ground). Parameters are checked by the :class:`.CuboidLeafCloudParams` class, which allows for many parameter combinations. The produced leaf cloud uniformly covers the :math:`(x, y, z) \\in \\left[ -\\dfrac{l_h}{2}, + \\dfrac{l_h}{2} \\right] \\times \\left[ -\\dfrac{l_h}{2}, + \\dfrac{l_h}{2} \\right] \\times [0, l_v]` region. Leaf orientation is controlled by the ``mu`` and ``nu`` parameters of an approximated inverse beta distribution :cite:`Ross1991MonteCarloMethods`. Finally, extra parameters control the random number generator and a basic and conservative leaf collision detection algorithm. Parameters ---------- seed : int Seed for the random number generator. avoid_overlap : bool If ``True``, generate leaf positions with strict collision checks to avoid overlapping. n_attempts : int If ``avoid_overlap`` is ``True``, number of attempts made at placing a leaf without collision before giving up. Default: 1e5. **kwargs Keyword arguments interpreted by :class:`.CuboidLeafCloudParams`. Returns ------- .LeafCloud Generated leaf cloud. See Also -------- :class:`.CuboidLeafCloudParams` """ rng = np.random.default_rng(seed=seed) n_attempts = kwargs.pop("n_attempts", int(1e5)) params = CuboidLeafCloudParams(**kwargs) if avoid_overlap: leaf_positions = _leaf_cloud_positions_cuboid_avoid_overlap( params.n_leaves, params.l_horizontal, params.l_vertical, params.leaf_radius, n_attempts, rng, ) else: leaf_positions = _leaf_cloud_positions_cuboid( params.n_leaves, params.l_horizontal, params.l_vertical, rng ) leaf_orientations = _leaf_cloud_orientations( params.n_leaves, params.mu, params.nu, rng ) leaf_radii = _leaf_cloud_radii(params.n_leaves, params.leaf_radius) # Create leaf cloud object return cls( id=params.id, leaf_positions=leaf_positions, leaf_orientations=leaf_orientations, leaf_radii=leaf_radii, leaf_reflectance=params.leaf_reflectance, leaf_transmittance=params.leaf_transmittance, )
[docs] @classmethod def sphere(cls, seed: int = 12345, **kwargs) -> LeafCloud: """ Generate a leaf cloud with spherical shape. Parameters are checked by the :class:`.SphereLeafCloudParams` class. The produced leaf cloud covers uniformly the :math:`r < \\mathtt{radius}` region. Leaf orientation is controlled by the ``mu`` and ``nu`` parameters of an approximated inverse beta distribution :cite:`Ross1991MonteCarloMethods`. An additional parameter controls the random number generator. Parameters ---------- seed : int Seed for the random number generator. **kwargs Keyword arguments interpreted by :class:`.SphereLeafCloudParams`. Returns ------- :class:`.LeafCloud` Generated leaf cloud. See Also -------- :class:`.SphereLeafCloudParams` """ rng = np.random.default_rng(seed=seed) params = SphereLeafCloudParams(**kwargs) leaf_positions = _leaf_cloud_positions_ellipsoid( params.n_leaves, rng, params.radius, params.radius, params.radius ) leaf_orientations = _leaf_cloud_orientations( params.n_leaves, params.mu, params.nu, rng ) leaf_radii = _leaf_cloud_radii(params.n_leaves, params.leaf_radius) # Create leaf cloud object return cls( id=params.id, leaf_positions=leaf_positions, leaf_orientations=leaf_orientations, leaf_radii=leaf_radii, leaf_reflectance=params.leaf_reflectance, leaf_transmittance=params.leaf_transmittance, )
[docs] @classmethod def ellipsoid(cls, seed: int = 12345, **kwargs) -> LeafCloud: """ Generate a leaf cloud with ellipsoid shape. Parameters are checked by the :class:`.EllipsoidLeafCloudParams` class. The produced leaf cloud covers uniformly the volume enclosed by :math:`\\frac{x^2}{a^2} + \\frac{y^2}{b^2} + \\frac{z^2}{c^2}= 1` . Leaf orientation is controlled by the ``mu`` and ``nu`` parameters of an approximated inverse beta distribution :cite:`Ross1991MonteCarloMethods`. An additional parameter controls the random number generator. Parameters ---------- seed : int Seed for the random number generator. **kwargs Keyword arguments interpreted by :class:`.EllipsoidLeafCloudParams`. Returns ------- :class:`.LeafCloud` Generated leaf cloud. See Also -------- :class:`.EllipsoidLeafCloudParams` """ rng = np.random.default_rng(seed=seed) params = EllipsoidLeafCloudParams(**kwargs) leaf_positions = _leaf_cloud_positions_ellipsoid( params.n_leaves, rng, params.a, params.b, params.c ) leaf_orientations = _leaf_cloud_orientations( params.n_leaves, params.mu, params.nu, rng ) leaf_radii = _leaf_cloud_radii(params.n_leaves, params.leaf_radius) # Create leaf cloud object return cls( id=params.id, leaf_positions=leaf_positions, leaf_orientations=leaf_orientations, leaf_radii=leaf_radii, leaf_reflectance=params.leaf_reflectance, leaf_transmittance=params.leaf_transmittance, )
[docs] @classmethod def cylinder(cls, seed: int = 12345, **kwargs) -> LeafCloud: """ Generate a leaf cloud with a cylindrical shape (vertical orientation). Parameters are checked by the :class:`.CylinderLeafCloudParams` class. The produced leaf cloud covers uniformly the :math:`r < \\mathtt{radius}, z \\in [0, l_v]` region. Leaf orientation is controlled by the ``mu`` and ``nu`` parameters of an approximated inverse beta distribution :cite:`Ross1991MonteCarloMethods`. An additional parameter controls the random number generator. Parameters ---------- seed : int Seed for the random number generator. **kwargs Keyword arguments interpreted by :class:`.CylinderLeafCloudParams`. Returns ------- :class:`.LeafCloud` Generated leaf cloud. See Also -------- :class:`.CylinderLeafCloudParams` """ rng = np.random.default_rng(seed=seed) params = CylinderLeafCloudParams(**kwargs) leaf_positions = _leaf_cloud_positions_cylinder( params.n_leaves, params.radius, params.l_vertical, rng ) leaf_orientations = _leaf_cloud_orientations( params.n_leaves, params.mu, params.nu, rng ) leaf_radii = _leaf_cloud_radii(params.n_leaves, params.leaf_radius) # Create leaf cloud object return cls( id=params.id, leaf_positions=leaf_positions, leaf_orientations=leaf_orientations, leaf_radii=leaf_radii, leaf_reflectance=params.leaf_reflectance, leaf_transmittance=params.leaf_transmittance, )
[docs] @classmethod def cone(cls, seed: int = 12345, **kwargs) -> LeafCloud: """ Generate a leaf cloud with a right conical shape (vertical orientation). Parameters are checked by the :class:`.ConeLeafCloudParams` class. The produced leaf cloud covers uniformly the :math:`r < \\mathtt{radius} \\cdot \\left( 1 - \\frac{z}{l_v} \\right), z \\in [0, l_v]` region. Leaf orientation is controlled by the ``mu`` and ``nu`` parameters of an approximated inverse beta distribution :cite:`Ross1991MonteCarloMethods`. An additional parameter controls the random number generator. Parameters ---------- seed : int Seed for the random number generator. **kwargs Keyword arguments interpreted by :class:`.ConeLeafCloudParams`. Returns ------- :class:`.LeafCloud` Generated leaf cloud. See Also -------- :class:`.ConeLeafCloudParams` """ rng = np.random.default_rng(seed=seed) params = ConeLeafCloudParams(**kwargs) leaf_positions = _leaf_cloud_positions_cone( params.n_leaves, params.radius, params.l_vertical, rng ) leaf_orientations = _leaf_cloud_orientations( params.n_leaves, params.mu, params.nu, rng ) leaf_radii = _leaf_cloud_radii(params.n_leaves, params.leaf_radius) # Create leaf cloud object return cls( id=params.id, leaf_positions=leaf_positions, leaf_orientations=leaf_orientations, leaf_radii=leaf_radii, leaf_reflectance=params.leaf_reflectance, leaf_transmittance=params.leaf_transmittance, )
[docs] @classmethod def from_file( cls, filename, leaf_transmittance: float | Spectrum = 0.5, leaf_reflectance: float | Spectrum = 0.5, id: str = "leaf_cloud", ) -> LeafCloud: """ Construct a :class:`.LeafCloud` from a text file specifying the leaf positions and orientations. .. admonition:: File format Each line defines a single leaf with the following 7 numerical parameters separated by one or more spaces: * leaf radius; * leaf center (x, y and z coordinates); * leaf orientation (x, y and z of normal vector). .. important:: All quantities are assumed to be given in metre. Parameters ---------- filename : path-like Path to the text file specifying the leaves in the leaf cloud. Can be absolute or relative. leaf_reflectance : :class:`.Spectrum` or float Reflectance spectrum of the leaves in the cloud. Must be a reflectance spectrum (dimensionless). Default: 0.5. leaf_transmittance : :class:`.Spectrum` of float Transmittance spectrum of the leaves in the cloud. Must be a transmittance spectrum (dimensionless). Default: 0.5. id : str ID of the created :class:`.LeafCloud` instance. Returns ------- :class:`.LeafCloud`: Generated leaf cloud. Raises ------ Raises ------ FileNotFoundError If ``filename`` does not point to an existing file. """ if not os.path.isfile(filename): raise FileNotFoundError(f"no file at {filename} found.") radii_ = [] positions_ = [] orientations_ = [] with open(os.path.abspath(filename), "r") as definition_file: for i, line in enumerate(definition_file): values = [float(x) for x in line.split()] radii_.append(values[0]) positions_.append(values[1:4]) orientations_.append(values[4:7]) radii = np.array(radii_) * ureg.m positions = np.array(positions_) * ureg.m orientations = np.array(orientations_) return cls( id=id, leaf_positions=positions, leaf_orientations=orientations, leaf_radii=radii, leaf_reflectance=leaf_reflectance, leaf_transmittance=leaf_transmittance, )
# -------------------------------------------------------------------------- # Kernel dictionary generation # -------------------------------------------------------------------------- @property def bsdf_id(self) -> str: return f"bsdf_{self.id}" @property def _template_bsdfs(self) -> dict: objects = { "reflectance": traverse(self.leaf_reflectance)[0].data, "transmittance": traverse(self.leaf_transmittance)[0].data, } result = {f"{self.bsdf_id}.type": "bilambertian"} for obj_key, obj_template in objects.items(): for key, param in obj_template.items(): result[f"{self.bsdf_id}.{obj_key}.{key}"] = param return result @property def _template_shapes(self) -> dict: length_units = uck.get("length") result = {} bsdf_dict = {"type": "ref", "id": self.bsdf_id} for i_leaf, (position, normal, radius) in enumerate( zip( self.leaf_positions.m_as(length_units), self.leaf_orientations, self.leaf_radii.m_as(length_units), ) ): _, up = mi.coordinate_system(normal) to_world = mi.ScalarTransform4f.look_at( origin=position, target=position + normal, up=up ) @ mi.ScalarTransform4f.scale(radius) result[f"{self.id}_leaf_{i_leaf}"] = { "type": "disk", "bsdf": bsdf_dict, "to_world": to_world, } return result @property def _params_bsdfs(self) -> dict: objects = { "reflectance": traverse(self.leaf_reflectance)[1].data, "transmittance": traverse(self.leaf_transmittance)[1].data, } result = {} for obj_key, obj_template in objects.items(): for key, param in obj_template.items(): # If no lookup strategy is set, we must add one if isinstance(param, UpdateParameter) and param.lookup_strategy is None: param = attrs.evolve( param, lookup_strategy=TypeIdLookupStrategy( mi.BSDF, self.bsdf_id, parameter_relpath=f"{obj_key}.{key}", ), ) result[f"{self.bsdf_id}.{obj_key}.{key}"] = param return result @property def _params_shapes(self) -> dict: return {} # -------------------------------------------------------------------------- # Other methods # --------------------------------------------------------------------------
[docs] def translated(self, xyz: pint.Quantity) -> LeafCloud: """ Return a copy of self translated by the vector ``xyz``. Parameters ---------- xyz : :class:`pint.Quantity` A 3-vector or a (N, 3)-array by which leaves will be translated. If (N, 3) variant is used, the array shape must match that of ``leaf_positions``. Returns ------- :class:`LeafCloud` Translated copy of self. Raises ------ ValueError Sizes of ``xyz`` and ``self.leaf_positions`` are incompatible. """ if xyz.ndim <= 1: xyz = xyz.reshape((1, 3)) elif xyz.shape != self.leaf_positions.shape: raise ValueError( f"shapes xyz {xyz.shape} and self.leaf_positions " f"{self.leaf_positions.shape} do not match" ) return attrs.evolve(self, leaf_positions=self.leaf_positions + xyz)