# Copyright 2021 United States Government as represented by the Administrator of the National Aeronautics and Space
# Administration. No copyright is claimed in the United States under Title 17, U.S. Code. All Other Rights Reserved.
"""
This module defines the illumination functions that are used to convert ray traced geometry into intensity values when
rendering.
In GIANT, an illumination model is typically a callable that converts an array of ray traced geometry with dtype
:attr:`.ILLUM_DTYPE` (typically generated by :meth:`.Scene.get_illumination_inputs`) into an array of floats giving the
intensity that would be expected to be returned along that ray. You can then add these values to a 2D image to create
a rendered image. Note that this does not take into account camera electronics or actual photon counts or anything
like that, since this is typically not important in OpNav. You could conceivably roll your own solution using the
provided components to do something like that if you so desired.
The typical workflow goes as follows:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from giant.ray_tracer.illumination import McEwenIllumination
>>> from giant.ray_tracer.scene import Scene, SceneObject
>>> from giant.ray_tracer.shapes import Ellipsoid, Point
>>> from giant.ray_tracer.rays import compute_rays
>>> from giant.camera_models import PinholeModel
>>> # define the camera model
>>> model = PinholeModel(focal_length=7.2, kx=1/2.2e-3, ky=1/2.2e-3, n_rows=1024, n_cols=1024, px=511.5, py=511.5)
>>> # set up the target that we want to render
>>> target = SceneObject(Ellipsoid(principal_axes=[0.25, 0.2, 0.25]))
>>> # move the target along the camera frame z axis
>>> target.translate([0, 0, 2.5])
>>> # create the sun
>>> sun = SceneObject(Point([0, 0, 0]))
>>> # put the sun behind and to the right of the camera
>>> sun.translate([800, 0, -400])
>>> # create the scene
>>> scene = Scene(target_objs=[target], light_obj=sun)
>>> # get the rays to trace through the scene (for the entire image)
>>> rays, pix = compute_rays(model, (0, model.n_rows-1), (0, model.n_cols-1))
>>> # get the illumination inputs using the scene and the rays
>>> illum_inputs = scene.get_illumination_inputs(rays)
>>> # use the illumination function to convert the illumination inputs into intensity values
>>> intensity = McEwenIllumination()(illum_inputs)
>>> # create the image array
>>> image = np.zeros((model.n_rows, model.n_cols), dtype=np.float64)
>>> # set up the subscripts where we need to add the intensity values to the image
>>> subs = np.round(pix).astype(int)
>>> # add the intensity values to the image
>>> np.add.at(image, (subs[1], subs[0]), intensity.ravel())
>>> # show it
>>> plt.imshow(image, cmap='gray')
>>> plt.show()
The most commonly used illumination function for OpNav purposes is the :class:`.McEwenIllumination`, which is a linear
weighting of the Lommel-Seeliger and Lambertian functions. All of the illumination models in this module are classes,
which when initialized create callables that are then used for the conversion. A couple also provide other methods, as
discussed, which may be useful in some specific circumstances.
"""
from typing import Tuple
from abc import ABCMeta, abstractmethod
import numpy as np
ILLUM_DTYPE: np.dtype = np.dtype([("incidence", np.float64, (3,)),
("exidence", np.float64, (3,)),
("normal", np.float64, (3,)),
("albedo", np.float64), ("visible", bool)])
"""
The numpy datatype expected by the illumination functions in this module as input for conversion to intensity values.
For an overview of how structured data types work in numpy, refer to https://numpy.org/doc/stable/user/basics.rec.html
The following table describes the purpose of each field.
================ ================ ======================================================================================
Field Type Description
================ ================ ======================================================================================
incidence 3 element double The unit vector from the light source to the surface as a 3 element array
exidence 3 element double The unit vector from the surface to the camera as a 3 element array
normal 3 element double The unit vector from that is normal to the surface as a 3 element array
albedo double The albedo of the surface as a float
visible bool A boolean flag specifying if these components are visible. Something is typically
considered not visible if (a) nothing was struck by the ray the record corresponds to,
(b) the spot struck was in shadow, or (c) the spot struct had the normal vector
pointing in the wrong direction.
================ ================ ======================================================================================
"""
_EYE3: np.ndarray = np.eye(3, dtype=np.float64)
"""
Predefined 3x3 double precision identity matrix for optimization purposes
"""
_HXHYARR = np.array([[-1, 0], [0, -1], [0, 0]], dtype=np.float64)
"""
An array used in computing the slope jacobian for the :class:`.McEwenIllumination`
"""
[docs]class IlluminationModel(metaclass=ABCMeta):
"""
This abstract base class specifies the minimum interface expected of an illumination model in GIANT.
Currently this just involves the ability to call the instance given an array as input and getting an array as output
but in the future it may specify other required interfaces like Jacobian matrices.
"""
@abstractmethod
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
IlluminationModel subclass instances should be callable to generate illumination values given an input of the
normal, incidence, exidence, albedo values, and visible flag in the form of a structured numpy array
(:attr:`.ILLUM_DTYPE`).
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
pass
[docs]class LambertianIllumination(IlluminationModel):
r"""
This basic illumination model computes the intensity values as simply the cosine of the incidence angle times the
albedo.
Mathematically this is given by
.. math::
I = -\alpha_0\alpha\mathbf{n}^T\mathbf{i}
where :math:`\alpha_0` is the global albedo value, :math:`\alpha` is the local albedo, :math:`\mathbf{n}` is the
unit normal vector, :math:`\mathbf{i}` is the unit incidence vector, and :math:`I` is the intensity value.
This class makes use of a global albedo value stored in :attr:`global_albedo` which can be used to scale all outputs
from this model. For most cases of OpNav, however, this can be ignored.
Anywhere that is not visible (either because the ``visible`` flag was set to False, or the cosine of the incidence
angle is less than 0) is set to 0 in the return.
"""
def __init__(self, global_albedo: float = 1.):
"""
:param global_albedo: The global albedo used to scale all outputs.
"""
self.global_albedo: float = global_albedo
"""
The global albedo used to scale all output
"""
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
Computes the intensity using the Lambertian law.
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
# compute the cosine of the incidence angle
cos_inc_angle = -(illum_inputs["normal"] * illum_inputs["incidence"]).sum(axis=-1)
# initialize the illumination values array
illum_values = np.zeros(illum_inputs.shape, dtype=np.float64)
# check which observations are actually visible
visible_check = illum_inputs["visible"] & (cos_inc_angle >= 0)
# compute the illumination values for only the valid observations
illum_values[visible_check] = (self.global_albedo * illum_inputs[visible_check]["albedo"] *
cos_inc_angle[visible_check])
return illum_values
[docs]class LommelSeeligerIllumination(IlluminationModel):
r"""
This basic illumination model computes the intensity values as simply the cosine of the incidence angle divided by
the cosine of the incidence angle plus the exidence angle times the albedo.
Mathematically this is given by
.. math::
I = \alpha_0\alpha\frac{-\mathbf{n}^T\mathbf{i}}{-\mathbf{n}^T\mathbf{i}+\mathbf{n}^T\mathbf{e}}
where :math:`\alpha_0` is the global albedo value, :math:`\alpha` is the local albedo, :math:`\mathbf{n}` is the
unit normal vector, :math:`\mathbf{i}` is the unit incidence vector, :math:`\mathbf{e}` is the unit exidence vector,
and :math:`I` is the intensity value.
This class makes use of a global albedo value stored in :attr:`global_albedo` which can be used to scale all outputs
from this model. For most cases of OpNav, however, this can be ignored.
Anywhere that is not visible (either because the ``visible`` flag was set to False, or the cosine of the incidence
angle is less than 0) is set to 0 in the return.
"""
def __init__(self, global_albedo: float = 1.):
"""
:param global_albedo: The global albedo used to scale all outputs.
"""
self.global_albedo: float = global_albedo
"""
The global albedo used to scale all outputs
"""
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
Computes the intensity using the Lommel-Seeliger law.
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
# compute the cosine of the incidence and exidence angles
cos_inc_angle = -(illum_inputs["normal"] * illum_inputs["incidence"]).sum(axis=-1)
cos_emi_angle = (illum_inputs["normal"] * illum_inputs["exidence"]).sum(axis=-1)
# initialize the illumination output array
illum_values = np.zeros(illum_inputs.shape, dtype=np.float64)
# check which observations are actually valid
visible_check = illum_inputs["visible"] & (cos_inc_angle >= 0) & (cos_emi_angle >= 0)
# compute the illumination only for the valid points
illum_values[visible_check] = (self.global_albedo * illum_inputs[visible_check]["albedo"] *
cos_inc_angle[visible_check] /
(cos_inc_angle[visible_check] + cos_emi_angle[visible_check]))
return illum_values
[docs]class McEwenIllumination(IlluminationModel):
r"""
This illumination model computes the intensity values as the weighted sum between the Lommel-Seeliger and
Lambertian models, weighted using the phase angle.
Mathematically this is given by
.. math::
I = \alpha_0\alpha\left(-(1-\beta)\mathbf{n}^T\mathbf{i} +
2\beta\frac{-\mathbf{n}^T\mathbf{i}}{-\mathbf{n}^T\mathbf{i}+\mathbf{n}^T\mathbf{e}}\right)
where :math:`\alpha_0` is the global albedo value, :math:`\alpha` is the local albedo, :math:`\mathbf{n}` is the
unit normal vector, :math:`\mathbf{i}` is the unit incidence vector, :math:`\mathbf{e}` is the unit exidence vector,
:math:`\beta=e^{\frac{-a}{60}}`, :math:`a=\cos^{-1}{-\mathbf{i}^T\mathbf{e}}` is the phase angle in degrees,
and :math:`I` is the intensity value.
This class makes use of a global albedo value stored in :attr:`global_albedo` which can be used to scale all outputs
from this model. For most cases of OpNav, however, this can be ignored.
Anywhere that is not visible (either because the ``visible`` flag was set to False, or the cosine of the incidence
angle is less than 0) is set to 0 in the return.
This class also provides methods for both analytic and numeric Jacobian matrices, which can be using as part of
photoclinometry to estimate surfaces given intensity values.
"""
def __init__(self, global_albedo: float = 1.):
"""
:param global_albedo: The global albedo used to scale all outputs.
"""
self.global_albedo: float = global_albedo
"""
The global albedo used to scale all outputs
"""
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
Computes the intensity using the McEwen law.
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
# compute the cosine of the incidence and emission angles
cos_inc_angle = -(illum_inputs["normal"] * illum_inputs["incidence"]).sum(axis=-1)
cos_emi_angle = (illum_inputs["normal"] * illum_inputs["exidence"]).sum(axis=-1)
# compute the phase angle in degrees
phase_angle = np.arccos(-(illum_inputs["incidence"] * illum_inputs["exidence"]).sum(axis=-1)) * 180 / np.pi
# initialize the illumination array
illum_values = np.zeros(illum_inputs.shape, dtype=np.float64)
# determine which observations are actually valid
visible_check = illum_inputs["visible"] & (cos_inc_angle >= 0) & (cos_emi_angle >= 0)
# compute the weighting between the lommel seeliger and the lambertian
beta = np.exp(-phase_angle[visible_check] / 60.)
# compute the illumination values at only the valid observations
illum_values[visible_check] = (self.global_albedo * illum_inputs[visible_check]["albedo"] *
((1 - beta) * cos_inc_angle[visible_check] +
2 * beta * cos_inc_angle[visible_check] /
(cos_inc_angle[visible_check] + cos_emi_angle[visible_check])))
return illum_values
def _compute_dillum_dcosi(self, albedo: np.ndarray, beta: np.ndarray, cosi_dcosi_cose2: np.ndarray) -> np.ndarray:
r"""
Compute the change in the illumination with respect to a change in the cosine of the incidence angle
Mathematically this is given by
.. math::
\frac{\partial I}{\partial \cos(i)} \alpha_0\alpha(1-\beta+2\beta\frac{\cos{i}}{(\cos{i}+\cos{e})^2})
where :math:`\cos(i)` is the cosine of the incidence angle, :math:`\cos(e)` is the cosine of the exidence angle,
and all else is as defined before.
:param albedo: The albedo values as a numpy array of doubles
:param beta: The beta parameter as a numpy array of doubles
:param cosi_dcosi_cose2: The cosine of the incidence angle divided by the square of the sum of the cosine of the
incidence angle and the cosine of the exidence angle
:return: the change in the illumination given a change in the cosine of the incidence angle
"""
return self.global_albedo * albedo * (1 - beta + 2 * beta * cosi_dcosi_cose2)
def _compute_dillum_dcose(self, albedo: np.ndarray, beta: np.ndarray, cosi_dcosi_cose2: np.ndarray) -> np.ndarray:
r"""
Compute the change in the illumination with respect to a change in the cosine of the exidence angle
Mathematically this is given by
.. math::
\frac{\partial I}{\partial \cos(e)} = \alpha_0\alpha(-2\beta\frac{\cos{i}}{(\cos{i}+\cos{e})^2})
where :math:`\cos(i)` is the cosine of the incidence angle, :math:`\cos(e)` is the cosine of the exidence angle,
and all else is as defined before.
:param albedo: The albedo values as a numpy array of doubles
:param beta: The beta parameter as a numpy array of doubles
:param cosi_dcosi_cose2: The cosine of the incidence angle divided by the square of the sum of the cosine of the
incidence angle and the cosine of the exidence angle
:return: the change in the illumination given a change in the cosine of the exidence angle
"""
return -2 * beta * self.global_albedo * albedo * cosi_dcosi_cose2
@staticmethod
def _compute_dnhat_dn(normal: np.ndarray) -> np.ndarray:
r"""
Compute the change in the unit normal vector with respect to a change in the normal vector
Mathematically this is given by
.. math::
\frac{\partial\hat{\mathbf{n}}}{\partial\mathbf{n}} = \fac{\mathbf{n}\mathbf{n}^T-
\mathbf{n}^T\mathbf{n}\mathbf{I}}{(\mathbf{n}^T\mathbf{n})^{1.5}}
:param normal: the normal vector
:return: the change in the unit vector with respect to a change in the normal vector
"""
ntn = (normal ** 2).sum(axis=-1).reshape(-1, 1, 1)
ntn32 = ntn ** 1.5
otp = np.einsum('ij,jk->jik', normal.T, normal)
return (otp - _EYE3 * ntn) / ntn32
def _compute_dillum_dalbedo(self, cosi: np.ndarray, beta: np.ndarray, cosi_cose: np.ndarray) -> np.ndarray:
r"""
Compute the change in the illumination with respect to a change in the local albedo
Mathematically this is given by
.. math::
\frac{\partial I}{\partial \alpha} = \alpha_0\left(-(1-\beta)\mathbf{n}^T\mathbf{i} +
2\beta\frac{-\mathbf{n}^T\mathbf{i}}{-\mathbf{n}^T\mathbf{i}+\mathbf{n}^T\mathbf{e}}\right)
:param cosi: The cosine of the incidence angle
:param beta: the beta parameter
:param cosi_cose: the cosine of the incidence angle plus the cosine of the exidence angle
:return: the change int he illumination with respect to a change in the local albedo
"""
beta_cosi = beta*cosi
return self.global_albedo * (cosi - beta_cosi + 2 * beta_cosi / cosi_cose)
[docs] def compute_photoclinometry_jacobian(self, observations: np.ndarray, rotation_to_inertial: np.ndarray,
max_inc: float = 70 * np.pi / 180,
max_emi: float = 70 * np.pi / 180,
max_phase: float = 140 * np.pi / 180,
update: bool = False,
update_weight: float = 5e-3) -> Tuple[np.ndarray, np.ndarray]:
r"""
This computes the Jacobian matrix of the change in the illumination values given a change in the surface normal
(represented by a change in the surface slope) and a change in the local albedo values.
Mathematically the jacobian is
.. math::
\mathbf{J}=\frac{\partial I}{\partial\left[\begin{array} {ccc} h_x & h_y & \alpha \end{array}\right]} =
\left[\begin{array}{ccc} \frac{\partial I}{\partial h_x} & \frac{\partial I}{\partial h_y} &
\frac{\partial I}{\partial\alpha} \end{array}\right]
where
.. math::
\frac{\partial I}{\partial h_x} = \frac{\partial I}{\partial \hat{\mathbf{n}}}
\frac{\partial \hat{\mathbf{n}}}{\partial\mathbf{n}}\frac{\partial \mathbf{n}}{\partial h_x} \\
\frac{\partial I}{\partial h_y} = \frac{\partial I}{\partial \hat{\mathbf{n}}}
\frac{\partial \hat{\mathbf{n}}}{\partial\mathbf{n}}\frac{\partial \mathbf{n}}{\partial h_y} \\
\frac{\partial I}{\partial \alpha} = \alpha_0\left(-(1-\beta)\mathbf{n}^T\mathbf{i} +
2\beta\frac{-\mathbf{n}^T\mathbf{i}}{-\mathbf{n}^T\mathbf{i}+\mathbf{n}^T\mathbf{e}}\right) \\
\frac{\partial I}{\partial \hat{\mathbf{n}}} = \frac{\partial I}{\partial \cos(i)}
\frac{\partial \cos(i)}{\partial\hat{\mathbf{n}}} + \frac{\partial I}{\partial \cos(e)}
\frac{\partial \cos(e)}{\partial\hat{\mathbf{n}}}\\
\frac{\partial I}{\partial \cos(i)} = \alpha_0\alpha(-2\beta\frac{\cos{i}}{(\cos{i}+\cos{e})^2}) \\
\frac{\partial I}{\partial \cos(e)} = \alpha_0\alpha(-2\beta\frac{\cos{i}}{(\cos{i}+\cos{e})^2}) \\
\frac{\partial \cos(i)}{\partial \hat{\mathbf{n}}} = -\mathbf{i}^T \\
\frac{\partial \cos(e)}{\partial \hat{\mathbf{n}}} = \mathbf{e}^T \\
\frac{\partial\hat{\mathbf{n}}}{\partial\mathbf{n}} = \frac{\mathbf{n}\mathbf{n}^T-
\mathbf{n}^T\mathbf{n}\mathbf{I}}{(\mathbf{n}^T\mathbf{n})^{1.5}}
:math:`\cos(i)` is the cosine of the incidence angle, :math:`\cos(e)` is the cosine of the exidence angle,
and all else is as defined before.
:param observations: The observations as a numpy array with type :attr:`.ILLUM_DTYPE`.
:param rotation_to_inertial: The rotation that takes the frame the observations are expressed in into the the
local frame for the surface (usually the local east north up frame)
:param max_inc: the maximum incidence angle to consider valid in radians
:param max_emi: the maximum emission angle to consider valid in radians
:param max_phase: the maximum phase angle to consider valid in radians
:param update: A flag specifying whether to form an update Jacobian
:param update_weight: The weight of the prior values if doing an update
:return: the jacobian matrix as a n(+3)x3 array (where n is the number of observations) and a boolean array of
length n specifying which rows of the jacobian matrix are valid
"""
if update:
out_jac = np.zeros((observations.shape[0] + 3, 3), dtype=np.float64)
out_jac[[-3, -2, -1], [-3, -2, -1]] = np.sqrt(update_weight)
else:
out_jac = np.zeros((observations.shape[0], 3), dtype=np.float64)
# extract the vectors we need from the structured array for ease of use
inc = -observations['incidence']
exi = observations['exidence']
nor = observations['normal']
alb = observations['albedo']
# compute the cosine of the incidence and exidence angles
cos_inc_angle = (nor * inc).sum(axis=-1)
cos_emi_angle = (nor * exi).sum(axis=-1)
# compute some common terms for efficiency
cosi_cose = cos_emi_angle + cos_inc_angle
cosi_cose2 = cosi_cose**2
cosi_dcosi_cose2 = cos_inc_angle/cosi_cose2
# compute the phase angle
phase_angle = np.arccos((inc * exi).sum(axis=-1))
# compute the weighting between the lambertian and the lommel seeliger
beta = np.exp(-phase_angle / 1.0471975511965976)
# compute the change in the illumination with respect to a change in the cosine of the incidence angle
dillum_dcosi = self._compute_dillum_dcosi(alb, beta, cosi_dcosi_cose2)
# compute the change in the illumination with respect to a change in the cosine of the exidence angle
dillum_dcose = self._compute_dillum_dcose(alb, beta, cosi_dcosi_cose2)
# compute the change in the illumination with respect to a change in the albedo
dillum_da = self._compute_dillum_dalbedo(cos_inc_angle, beta, cosi_cose)
# the change in the cosine of the incidence angle with respect to a change in the normal vector
dcosi_dnhat = inc
# the change in the cosine of the exidence angle with respect to a change in the normal vector
dcose_dnhat = exi
# the change in the unit normal vector with respect to a change in the normal vector
dnhat_dn = self._compute_dnhat_dn(nor)
# the change in the normal vector with respect to a change in the slope
dn_dhxhy = rotation_to_inertial @ _HXHYARR
# the change in the illumination with respect to a change in the slope using chain rule
dillum_dhxhy = ((dillum_dcosi.reshape(-1, 1) * dcosi_dnhat +
dillum_dcose.reshape(-1, 1) * dcose_dnhat).reshape(-1, 1, 3) @ dnhat_dn @ dn_dhxhy).squeeze()
# initialize the boolean array for which observations are actually valid
valid = np.ones(out_jac.shape[0], dtype=bool)
if update:
# if we are doing an update step then we need to leave the last 3 rows alone
# check which observations are valid
valid[:-3] = (cos_inc_angle >= np.cos(max_inc)) & (cos_emi_angle >= np.cos(max_emi)) & \
(phase_angle <= max_phase) & (cos_inc_angle > 0) & (cos_emi_angle > 0)
# form the jacobian matrix
out_jac[:-3, :2] = dillum_dhxhy.reshape(-1, 2)
out_jac[:-3, 2] = dillum_da.ravel()
else:
# check which observations are valid
valid[:] = (cos_inc_angle >= np.cos(max_inc)) & (cos_emi_angle >= np.cos(max_emi)) & \
(phase_angle <= max_phase) & (cos_inc_angle > 0) & (cos_emi_angle > 0)
# form the jacobian matrix
out_jac[:, :2] = dillum_dhxhy.reshape(-1, 2)
out_jac[:, 2] = dillum_da.ravel()
return out_jac, valid
[docs] def numeric_derivative(self, observations: np.ndarray, rotation_to_inertial: np.ndarray, delta: float = 1e-6,
max_inc: float = 70 * np.pi / 180, max_emi: float = 70 * np.pi / 180,
max_phase: float = 140 * np.pi / 180) -> Tuple[np.ndarray, np.ndarray]:
r"""
This computes the Jacobian matrix of the change in the illumination values given a change in the surface normal
(represented by a change in the surface slope) and a change in the local albedo values use numeric finite
differencing.
Mathematically the jacobian is
.. math::
\mathbf{J}=\frac{\partial I}{\partial\left[\begin{array} {ccc} h_x & h_y & \alpha \end{array}\right]} =
\left[\begin{array}{ccc} \frac{\partial I}{\partial h_x} & \frac{\partial I}{\partial h_y} &
\frac{\partial I}{\partial\alpha} \end{array}\right]
Here we compute this using numeric central finite differencing.
This is primarily used for verifying the analytic jacobian given by :meth:`.compute_photoclinometry_jacobian`
which should be preferred in most cases due to its speed/efficiency.
:param observations: The observations as a numpy array with type :attr:`.ILLUM_DTYPE`.
:param rotation_to_inertial: The rotation that takes the frame the observations are expressed in into the the
local frame for the surface (usually the local east north up frame)
:param delta: The size of the perturbation to use in the finite differencing
:param max_inc: the maximum incidence angle to consider valid in radians
:param max_emi: the maximum emission angle to consider valid in radians
:param max_phase: the maximum phase angle to consider valid in radians
:return: the jacobian matrix as a n(+3)x3 array (where n is the number of observations) and a boolean array of
length n specifying which rows of the jacobian matrix are valid
"""
# numerically differentiate the BRDF
# compute the cosine of the incidence and exidence angles at the nominal geometry
cos_inc_angle = -(observations["normal"] * observations["incidence"]).sum(axis=-1)
cos_emi_angle = (observations["normal"] * observations["exidence"]).sum(axis=-1)
# compute the phase angle at the nominal geometry in degrees
phase_angle = np.arccos(-(observations["incidence"] * observations["exidence"]).sum(axis=-1)) * 180 / np.pi
# compute the normal vector in the inertial frame
normal = (rotation_to_inertial.T @ observations["normal"].T).T
normal /= normal[:, 2, np.newaxis]
# do the forward differences
# make a copy of the observations table
observations_pert = observations.copy()
# perturb the albedo channel
observations_pert["albedo"] += delta
# compute the values for the perturbed viewing conditions
illum_alb_f = self(observations_pert)
# make a copy of the observations table
observations_pert = observations.copy()
# perturb the slope in the x direction (the normal vector)
normal_inertial_dx_pert = (rotation_to_inertial @ (normal + [delta, 0, 0]).T).T
normal_inertial_dx_pert /= np.linalg.norm(normal_inertial_dx_pert, axis=-1, keepdims=True)
observations_pert["normal"] = normal_inertial_dx_pert
# compute the values for the perturbed viewing conditions
illum_dx_f = self(observations_pert)
# make a copy of the observations table
observations_pert = observations.copy()
# perturb the slope in the y direction (the normal vector)
normal_inertial_dy_pert = (rotation_to_inertial @ (normal + [0, delta, 0]).T).T
normal_inertial_dy_pert /= np.linalg.norm(normal_inertial_dy_pert, axis=-1, keepdims=True)
observations_pert["normal"] = normal_inertial_dy_pert
# compute the values for the perturbed viewing conditions
illum_dy_f = self(observations_pert)
# do the backward differences
# make a copy of the observations table
observations_pert = observations.copy()
# perturb the albedo channel
observations_pert["albedo"] -= delta
# compute the values for the perturbed viewing conditions
illum_alb_b = self(observations_pert)
# compute the values for the perturbed viewing conditions
observations_pert = observations.copy()
# perturb the slope in the x direction (the normal vector)
normal_inertial_dx_pert = (rotation_to_inertial @ (normal - [delta, 0, 0]).T).T
normal_inertial_dx_pert /= np.linalg.norm(normal_inertial_dx_pert, axis=-1, keepdims=True)
observations_pert["normal"] = normal_inertial_dx_pert
# compute the values for the perturbed viewing conditions
illum_dx_b = self(observations_pert)
# compute the values for the perturbed viewing conditions
observations_pert = observations.copy()
# perturb the slope in the y direction (the normal vector)
normal_inertial_dy_pert = (rotation_to_inertial @ (normal - [0, delta, 0]).T).T
normal_inertial_dy_pert /= np.linalg.norm(normal_inertial_dy_pert, axis=-1, keepdims=True)
observations_pert["normal"] = normal_inertial_dy_pert
# compute the values for the perturbed viewing conditions
illum_dy_b = self(observations_pert)
# determine which observations are actually valid
valid = (cos_inc_angle >= np.cos(max_inc)) & (cos_emi_angle >= np.cos(max_emi)) & \
(phase_angle <= max_phase * 180 / np.pi) & (cos_inc_angle > 0) & (cos_emi_angle > 0)
# construct the finite differenced jacobian
return np.array([(illum_dx_f - illum_dx_b) / (2 * delta),
(illum_dy_f - illum_dy_b) / (2 * delta),
(illum_alb_f - illum_alb_b) / (2 * delta)]).T, valid
[docs]class GaskellIllumination(IlluminationModel):
r"""
This illumination model computes the intensity values as the weighted sum between the Lommel-Seeliger and
Lambertian models, weighted using the phase angle.
Mathematically this is given by
.. math::
I = \alpha_0\alpha\left(-(1-\beta)\mathbf{n}^T\mathbf{i} +
\beta\frac{-\mathbf{n}^T\mathbf{i}}{-\mathbf{n}^T\mathbf{i}+\mathbf{n}^T\mathbf{e}}\right)
where :math:`\alpha_0` is the global albedo value, :math:`\alpha` is the local albedo, :math:`\mathbf{n}` is the
unit normal vector, :math:`\mathbf{i}` is the unit incidence vector, :math:`\mathbf{e}` is the unit exidence vector,
:math:`\beta=e^{\frac{-a}{60}}`, :math:`a=\cos^{-1}{-\mathbf{i}^T\mathbf{e}}` is the phase angle in degrees,
and :math:`I` is the intensity value. Note that this differs from the :class:`.McEwenIllumination` model only in
factor of 2 on the Lommel-Seeliger term. This was a typo in the SPC source code and is only included for
consistency therein.
This class makes use of a global albedo value stored in :attr:`global_albedo` which can be used to scale all outputs
from this model. For most cases of OpNav, however, this can be ignored.
Anywhere that is not visible (either because the ``visible`` flag was set to False, or the cosine of the incidence
angle is less than 0) is set to 0 in the return.
"""
def __init__(self, global_albedo: float = 1.):
"""
:param global_albedo: The global albedo used to scale all outputs.
"""
self.global_albedo: float = global_albedo
"""
The global albedo used to scale all outputs
"""
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
Computes the intensity using the Gaskell law.
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
# compute the cosine of the incidence and exidence angles
cos_inc_angle = -(illum_inputs["normal"] * illum_inputs["incidence"]).sum(axis=-1)
cos_emi_angle = (illum_inputs["normal"] * illum_inputs["exidence"]).sum(axis=-1)
# compute the phase angle in degrees
phase_angle = np.arccos(-(illum_inputs["incidence"] * illum_inputs["exidence"]).sum(axis=-1)) * 180 / np.pi
# initialize the output illumination array
illum_values = np.zeros(illum_inputs.shape, dtype=np.float64)
# check which observations are valid
visible_check = illum_inputs["visible"] & (cos_inc_angle >= 0) & (cos_emi_angle >= 0)
# compute the weighting between the lommel seeliger and lambertian for the valid observations
beta = np.exp(-phase_angle[visible_check] / 60.)
# compute the illumination values for the valid observations
illum_values[visible_check] = (self.global_albedo * illum_inputs[visible_check]["albedo"] *
((1 - beta) * cos_inc_angle[visible_check] +
beta * cos_inc_angle[visible_check] /
(cos_inc_angle[visible_check] + cos_emi_angle[visible_check])))
return illum_values
[docs]class AshikhminShirleyDiffuseIllumination(IlluminationModel):
r"""
This illumination model computes the intensity values according to the Ashikhmin Shirley diffuse law
Mathematically this is given by
.. math::
I = -\mathbf{n}^T\mathbf{i}(1-(1+\frac{\mathbf{n}^T\mathbf{i}}{2})^5)(1-(1-\frac{\mathbf{n}^T\mathbf{e}}{2})^5)
where :math:`\mathbf{n}` is the
unit normal vector, :math:`\mathbf{i}` is the unit incidence vector, :math:`\mathbf{e}` is the unit exidence vector,
and :math:`I` is the intensity value.
Anywhere that is not visible (either because the ``visible`` flag was set to False, or the cosine of the incidence
angle is less than 0) is set to 0 in the return.
This illumination model has not been thoroughly tested for use with natural bodies.
"""
def __call__(self, illum_inputs: np.ndarray) -> np.ndarray:
"""
Computes the intensity using the Ashikhmin-Shirley Diffuse law.
:param illum_inputs: The illumination inputs as a structured numpy array
:return: The intensity values for the input array as a numpy double precision array
"""
# compute the cosine of the incidence and exidence angles
cos_inc_angle = -(illum_inputs["normal"] * illum_inputs["incidence"]).sum(axis=-1)
cos_emi_angle = (illum_inputs["normal"] * illum_inputs["exidence"]).sum(axis=-1)
# determine which observations are actually valid
visible_check = illum_inputs["visible"] & (cos_inc_angle >= 0) & (cos_emi_angle >= 0)
# initialize the output illumination array
illum_values = np.zeros(illum_inputs.shape, dtype=np.float64)
# compute the illumination using only the valid observations
illum_values[visible_check] = (cos_inc_angle[visible_check] *
((1 - np.power(1 - cos_inc_angle[visible_check] / 2, 5)) *
(1 - np.power(1 - cos_emi_angle[visible_check] / 2, 5))))
return illum_values
