Understanding Radio Materials
Radio materials define how objects scatter incident radio waves.They implement all necessary components to simulate the interactionbetween radio waves and objects composed of specific materials.
Modifying Parameters of Radio Materials
To show how to modify the parameters of radio materials, we start by loadinga scene that consists only of a single reflector.We also instantiate a transmitter and a receiver, each equipped with a singleantenna.
scene=load_scene(rt.scene.simple_reflector,merge_shapes=False)scene.add(Transmitter("tx",position=(-2.,0.,1)))scene.add(Receiver("rx",position=(2.,0.,1)))scene.tx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")scene.rx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")
We then change the radio material of the reflector to aRadioMaterial instance.We set the relative permittivity and conductivity of the material to specific values.
# Instantiate the radio materialmy_mat=RadioMaterial("my-mat",thickness=0.1,relative_permittivity=5.,conductivity=1)# Assign the radio material to the reflectorscene.objects["reflector"].radio_material=my_mat# To avoid confusion, discard the radio material initially loaded with the scenescene.remove("reflector-mat")# Print the material to visualize its parametersprint(scene.radio_materials)
{'my-mat':RadioMaterialeta_r=5.000sigma=1.000thickness=0.100scattering_coefficient=0.000xpd_coefficient=0.000}
We can now compute paths and print their gains for different values of the conductivity.
# Instantiate the path solversolver=PathSolver()conductivities=[1.,100.,10000.,1000000.]forsigmainconductivities:my_mat.conductivity=sigmapaths=solver(scene)# Paths coefficienta,tau=paths.a,paths.taua=a[0].numpy()+1j*a[1].numpy()a=np.squeeze(a,(0,1,2,3))tau=tau.numpy()tau=np.squeeze(tau,(0,1))print("Conductivity [S/m]:",sigma)fora_,tau_inzip(a,tau):print("\t Delay [ns]: ",tau_*1e9," Gain [dB]",10.*np.log10(np.square(np.abs(a_))))
Conductivity[S/m]:1.0Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-69.881165Conductivity[S/m]:100.0Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-57.55436Conductivity[S/m]:10000.0Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-56.460648Conductivity[S/m]:1000000.0Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-56.351566
The two paths correspond to the line-of-sight and reflected paths, in that order.We can see how the reflected path gain increases as the conductivity of the reflectoris set to higher values.
Calibrating Material Parameters Through Gradient Descent
We consider a simple example in which we aim to retrieve the conductivity of thea radio material through gradient descent.
To that aim, we define a second scene with the same setup, except that theconductivity is initialized to an arbitrary value and is trainable.We will then retrieve the value of the conductivity through gradient descent onthe normalized absolute error between the path gain computed for the previouslydefined scene (reference scene) and the trainable one.
trainable_scene=load_scene(rt.scene.simple_reflector,merge_shapes=False)trainable_scene.add(Transmitter("tx",position=(-2.,0.,1)))trainable_scene.add(Receiver("rx",position=(2.,0.,1)))trainable_scene.tx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")trainable_scene.rx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")
In order to assign trainable variables to the reflector material properties, weneed to first instantiate an optimizer.The optimizer is used to instantiate trainable variables.
# Adam optimizerlearning_rate=4e-2opt=mi.ad.Adam(lr=learning_rate)
We then replace the material reflector by a newly instantiated material with atrainable conductivity property.To ensure that the conductivity remains non-negative, the trainable variable isdefined as its logarithm (logit), and the actual conductivity is computed fromthe logit using a numerically stablesigmoid() function.
deflogit_2_conductivity(logit):max_conductivity=100.returnsigmoid(logit)*max_conductivity# Trainable variable is the logit of the conductivity# It is initialized to an arbitrary valueopt["logit_conductivity"]=mi.Float(-5.)# Instantiate the radio material# The conductivity is initialized using the trainable variabletrainable_mat=RadioMaterial("my-trainable-mat",thickness=0.1,relative_permittivity=5.,conductivity=logit_2_conductivity(opt["logit_conductivity"]))# Use the trainable variable# Assign the radio material to the reflectortrainable_scene.objects["reflector"].radio_material=trainable_mat# To avoid confusion, discard the radio material initially loaded with the scenetrainable_scene.remove("reflector-mat")# Print the material to visualize its parametersprint(trainable_scene.radio_materials)
{'my-trainable-mat':RadioMaterialeta_r=5.000sigma=0.669thickness=0.100scattering_coefficient=0.000xpd_coefficient=0.000}
Backpropagation through adrjit loop can currently only be done when theevaluated mode is used for the computation of the electric field by the solver.
solver=PathSolver()# Switch the computation of field loop to "evaluated" mode to# enable gradient backpropagation through the loopsolver.loop_mode="evaluated"
We can now run the optimization that performs gradient descent on the normalizedabsolute error between the gain obtained using the trainable scene and the oneobtained using the reference scene instantiated at the beginning of this guide.
deftotal_gain(paths):a_real,a_imag=paths.ag=dr.sum(dr.square(a_real)+dr.square(a_imag))returng# Normalized absolute error# The `dr.detach()` function stops gradient from# propagating through its inputdefnae(x,y):returndr.abs(x-y)*dr.detach(dr.rcp(y))num_iterations=150# Ground-truthscene.radio_materials["my-mat"].conductivity=20.ref_paths=solver(scene)ref_gain=dr.detach(total_gain(ref_paths))# Record the conductivity valueconductivity=[]# Optimization loopfor_inrange(num_iterations):# Run simulationpaths=solver(trainable_scene)# Compute loss on total gain and the gradientsgain=total_gain(paths)loss=nae(gain,ref_gain)dr.backward(loss)# Optimizer stepopt.step()updated_conductivity=logit_2_conductivity(opt["logit_conductivity"])conductivity.append(updated_conductivity.numpy())trainable_mat.conductivity=updated_conductivityplt.figure()plt.grid(True)plt.plot(np.arange(1,num_iterations+1),conductivity,label="Calibrated")plt.hlines([scene.get("my-mat").conductivity],1,num_iterations+1,color="k",label="Ground-truth")plt.xlabel("Iteration")plt.ylabel("Conductivity [S/m]")plt.legend();
Custom Radio Materials
Compared to what was done in the previous section, we will now implement ascattering model by defining a new class that inherits from theRadioMaterialBase abstract class, which allows us to freelydefine how a material scatters an incident wave.
We will start by detailing how Jones vectors and matrices are represented inSionna RT. It is highly recommended to first read thePrimer on Electromagnetics to understand the basics ofradio wave propagation.
Representation of Jones vector and Matrices
As detailed in thePrimer on Electromagnetics, a wave phasoris typically represented by a Jones vector\(\mathbf{E} \in \mathbb{C}^2\).However, asdrjit does not currently support complex-valued vectors andmatrices, Sionna RT uses the equivalent real-valued representation of two-dimensionalcomplex-valued vectors, which represents a Jones vector as a real-valued vectorwith 4 dimensions:
Similarly, a Jones matrix\(\mathbf{M}\), which is a complex-valued matrixof size\(2 \times 2\) that models an interaction with the environment, isrepresented by the equivalent\(4 \times 4\) real-valued matrix
In the rest of this guide, we will interchangeably use both representations towrite equations.However, all implementations use the real-valued representation.
Implicit Basis
A wave phasor\(\mathbf{E}\) is expressed by two arbitrary orthogonalpolarization directions S and P:
which are both orthogonal to the direction of propagation vector\(\hat{\mathbf{k}}\),i.e.,\(\hat{\mathbf{e}}_{s}^{\mathsf{T}} \hat{\mathbf{e}}_{p}=\hat{\mathbf{e}}_{s}^{\mathsf{T}} \hat{\mathbf{k}}=\hat{\mathbf{e}}_{p}^{\mathsf{T}} \hat{\mathbf{k}} = 0\).As\((\hat{\mathbf{e}}_{s}, \hat{\mathbf{e}}_{p}, \hat{\mathbf{k}})\) forms an orthonormal basis, only the knowledge of the S basis vector\(\hat{\mathbf{e}}_{s}\) (in addition tothe direction of propagation\(\hat{\mathbf{k}}\)) is required to define thereference basis in which the wave phasor is represented.
To avoid storing the S basis vector for every wave phasor, animplicit basis is used in Sionna RT.More precisely, a wave phasor that propagates in direction\(\hat{\mathbf{k}}\)is always assumed to use as S basis vector the unit vector computed from itsdirection of propagation by the utilityimplicit_basis_vector().This utility deterministically constructs a unit vector orthogonal to\(\hat{\mathbf{k}}\) to be used as S basis vector\(\hat{\mathbf{e}}_{s}\).
When implementing a custom radio material in Sionna RT, it is therefore requiredto compute the Jones matrix describing the interaction with the material assumingthat the Jones vector on which the matrix is applied describes the wave using theimplicit basis.Moreover, it is required that the result of applying this Jones matrix is aJones vector that also describes the scattered wave using the implicit basis.
The Local Interaction Basis
Computing the Jones matrix and direction of propagation of the scattered waveresulting from an interaction is facilitated in Sionna RT by defining a localcoordinate system specific to this interaction.
This local interaction coordinate system is defined such that the local Z axiscorresponds to the normal to the scatterer surface at the interaction point, andthe X and Y axes are tangent to the scatterer surface at the interaction point.
Using the prime (\('\)) notation to indicate that a vector is representedin the local coordinate system, we therefore have\(\mathbf{n'} = \begin{bmatrix}0\\0\\1\end{bmatrix}\).
Mandatory Subclass Methods
Implementing a custom radio material requires defining a class that inherits fromRadioMaterialBase and implements the following methods:
sample()
Samples an interaction type and the direction of propagation of the scattered wave(which typically depend on the sampled interaction type).This function must return, among others, the sample interaction type, directionof propagation of the scattered wave, and Jones matrix that models the interaction.
Sampling of the material is required as, when a ray modeling an incident waveintersects a scatterer, only a single ray can be spawn to model the scattered field.Indeed, spawning more than one ray per interaction (e.g., one refracted and onereflected ray, or multiple diffusely reflected rays) is computationallyprohibitive as it leads to an exponential scaling of the complexity with the path length.Therefore, a single ray corresponding to a single interaction type (e.g.,refracted or reflected) and a single direction of propagation are randomly sampled.Accurate modeling of the scattered field is achieved by having a large numberof incident rays interacting with the material resulting in independently sampledscattered rays that model well the scattered field.
eval()
Evaluates the Jones matrix for a given interaction type, direction of incidence,and direction of scattering. Compared tosample(),this method does not sample the material.
pdf()
Returns the probability that a given interaction type and direction of scatteringare sampled conditioned on a given direction of incidence.
traverse()
Traverses the attributes and objects of the material. This method is used torecord the material parameters, and especially the differentiable parameters.
to_string()
Returns a string describing the material. This is used to “print” the materialin a humanly readable way.
Implementation of a Simple Radio Material Model
For simplicity, we will start by implementing a scattering model that onlyreflects incident radio waves specularly, and such that the energy of the reflectedwave equals the one of the incident wave scaled by a parameter\(g \in (0,1)\).
Formally, this simple material model implements the following transfer equation.The direction of propagation of the reflected ray\(\hat{\mathbf{k}}_\text{r}\)is computed as:
where\(\hat{\mathbf{k}}_\text{i}\) is the direction of propagation of theincident wave and\(\hat{\mathbf{n}}\) the normal to the scatterer at theinteraction point.
The reflected wave phasor\(\mathbf{E}_\text{r}\) is represented as
where
and
is the incident wave phasor.In these equations,\((\hat{\mathbf{e}}_{\text{i},s}, \hat{\mathbf{e}}_{\text{i},p})\)and\((\hat{\mathbf{e}}_{\text{r},s}, \hat{\mathbf{e}}_{\text{r},p})\) arethe implicit basis for the incident and scattered wave, respectively,and\(\mathbf{W}\left(\hat{\mathbf{x}}_{2}, \hat{\mathbf{y}}_{2}, \hat{\mathbf{x}}_{1}, \hat{\mathbf{y}}_{1}\right)\)is the change-of-basis matrix from basis\((\hat{\mathbf{x}}_{1}, \hat{\mathbf{y}}_{1})\) tobasis\((\hat{\mathbf{x}}_{2}, \hat{\mathbf{y}}_{2})\) defined in thePrimer on Electromagnetics .
Let’s now implement this radio material model.Note that only thesample(),eval(),andpdf() methods are required to actuallyimplement the previous equations.
The direction of propagation of the incident wave is assumed to be provided tothe radio material represented in the local coordinate system.The direction of reflection(56) can then be computed by leveraging the fact that\(\hat{\mathbf{n}} = \hat{\mathbf{z}}\) as shown in the following code snippet, wherean arbitrary direction of incidence is defined.
# Arbitrary direction of incidence used in this example.# The z component of the incident direction of propagation can only be negativeki_prime=dr.normalize(mi.Vector3f(1.,2.,-1.))print("ki_prime",ki_prime)# By definitionn_prime=mi.Vector3f(0.,0.,1.)kr_prime=ki_prime-2.*dr.dot(n_prime,ki_prime)*n_primeprint("kr_prime",kr_prime)
ki_prime[[0.408248,0.816497,-0.408248]]kr_prime[[0.408248,0.816497,0.408248]]
Equivalently, one can use the Mitsuba utilitymi.reflect().Note that this utility uses the rendering convention for representing directionof propagation, in which the input vector pointsaway from the interactionpoint:
# Expects an input vector that points away from the interaction pointprint("kr_prime",mi.reflect(-ki_prime))
kr_prime[[0.408248,0.816497,0.408248]]
To compute the scattered field(57), the Jones matrix needs to beconstructed such that both the incident field and the scattered field (resultingfrom applying the Jones matrix) are represented in their respective implicitworld coordinate system.Therefore, the Jones matrix in(57), i.e.,
is such that the change-of-basis matrix\(\mathbf{W}\left(\hat{\mathbf{e}}_{\text{i},\perp}, \hat{\mathbf{e}}_{\text{i},\parallel}, \hat{\mathbf{e}}_{\text{i},s}, \hat{\mathbf{e}}_{\text{i},p}\right)\) is from the implicit world coordinate system for the direction of propagation\(\hat{\mathbf{k}}_i\) to the coordinate system where\(\hat{\mathbf{e}}_{\text{i},\perp}\)is the TE polarization and\(\hat{\mathbf{e}}_{\text{i},\parallel}\) is theTM polarization for the incident wave.Similarly,\(\mathbf{W}\left(\hat{\mathbf{e}}_{\text{r},s}, \hat{\mathbf{e}}_{\text{r},p}, \hat{\mathbf{e}}_{\text{r},\perp}, \hat{\mathbf{e}}_{\text{r},\parallel}\right)\) is the change-of-basis matrix from the coordinate system where\(\hat{\mathbf{e}}_{\text{r},\perp}\) is the TE polarization and\(\hat{\mathbf{e}}_{\text{r},\parallel}\) is the TM polarization for the scatteredwave to the implicit world coordinate system for the direction of propagation\(\hat{\mathbf{k}}_r\).
As Jones matrices with this structure are commonly used to model specular reflectionand refraction, Sionna RT provides thejones_matrix_to_world_implicit()utility to construct such matrices.In this case, we can utilize this function by setting bothc1 andc2 parameters to\(\sqrt{g}\).
The next code snippet provides an implementation of the previous model as a customradio material, with comments explaining every step.Consulting the API documentation ofRadioMaterialBase is essentialfor a detailed understanding of this code.It is also recommended to consult theAPI documentation of the mi.BSDF class,asRadioMaterialBase inherits from it.
classCustomRadioMaterial(RadioMaterialBase):# The __init__ method builds the radio material from:# - A unique `name` to identify the material instance in the scene# - The gain parameter `g`# - An optional `color` for displaying the material in the previewer and renderer# Providing these 3 parameters to __init__ is how an instance of this radio material# is built programmatically.## When loading a scene from an XML file, Mitsuba provides to __init__# only an `mi.Properties` object containing all the properties of the material# read from the XML scene file. Therefore, when a `props` object is provided,# the other parameters are ignored and should not be given.def__init__(self,name:str|None=None,g:float|mi.Float|None=None,color:Tuple[float,float,float]|None=None,props:mi.Properties|None=None):# If `props` is `None`, then one is built from the# other parametersifpropsisNone:props=mi.Properties("custom-radio-material")# Name of the radio materialprops.set_id(name)props["g"]=gifcolorisnotNone:props["color"]=mi.ScalarColor3f(color)# Read the gain from `props`g=0.0if"g"inprops:g=props["g"]delprops["g"]self._g=mi.Float(g)# The other parameters (`name`, `color`) are given to the# base class to complete the initialization of the materialsuper().__init__(props)defsample(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,sample1:mi.Float,sample2:mi.Point2f,active:bool|mi.Bool=True)->Tuple[mi.BSDFSample3f,mi.Spectrum]:# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Build the 3x3 change-of-basis matrix from the local basis to the world# basis.# `si.sh_frame` stores the three vectors that define the local interaction basis# in the world coordinate system.to_world=mi.Matrix3f(si.sh_frame.s,si.sh_frame.t,si.sh_frame.n).T# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# Compute the Jones matrix in the implicit world coordinate system# The `jones_matrix_to_world_implicit()` builds the Jones matrix with the# structure we need.sqrt_g=mi.Complex2f(dr.sqrt(self._g),0.)jones_mat=jones_matrix_to_world_implicit(c1=sqrt_g,c2=sqrt_g,to_world=to_world,k_in_local=ki_prime,k_out_local=kr_prime)## We now only need to prepare the outputs# Cast the Jones matrix to a `mi.Spectrum` to meet the requirements of# the BSDF interface of Mitsubajones_mat=mi.Spectrum(jones_mat)# Instantiate and set the BSDFSample objectbs=mi.BSDFSample3f()# Specifies the type of interaction that was computedbs.sampled_component=InteractionType.SPECULAR# Direction of the scattered wave in the world framebs.wo=to_world@kr_prime# The next field of `bs` stores the probability that the sampled# interaction type and direction of scattering are sampled conditioned# on the given direction of incidence.# As only one event and direction of scattering are possible with this model,# this probability is set to 1.bs.pdf=mi.Float(1.)# Not used but required to be setbs.sampled_type=mi.UInt32(+mi.BSDFFlags.DeltaReflection)bs.eta=1.0returnbs,jones_matdefeval(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,wo:mi.Vector3f,active:bool|mi.Bool=True)->mi.Spectrum:# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Build the 3x3 change-of-basis matrix from the local basis to the world# basis.# `si.sh_frame` stores the three vectors that define the local interaction basis# in the world coordinate system.to_world=mi.Matrix3f(si.sh_frame.s,si.sh_frame.t,si.sh_frame.n).T# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# Compute the Jones matrix in the implicit world coordinate system# The `jones_matrix_to_world_implicit()` builds the Jones matrix with the# structure we need.sqrt_g=mi.Complex2f(dr.sqrt(self._g),0.)jones_mat=jones_matrix_to_world_implicit(c1=sqrt_g,c2=sqrt_g,to_world=to_world,k_in_local=ki_prime,k_out_local=kr_prime)# This model only scatters energy in the direction of the specular reflection.# Any other direction provided by the user `wo` should therefore lead to no energy.is_valid=isclose(dr.dot(kr_prime,wo),mi.Float(1.))jones_mat=dr.select(is_valid,jones_mat,0.)# Cast the Jones matrix to a `mi.Spectrum` to meet the requirements of# the BSDF interface of Mitsubajones_mat=mi.Spectrum(jones_mat)returnjones_matdefpdf(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,wo:mi.Vector3f,active:bool|mi.Bool=True)->mi.Float:# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# As only one event and direction of scattering are possible with this model,# the probability is set to 1 for this direction and 0 for any other.is_valid=isclose(dr.dot(kr_prime,wo),mi.Float(1.))returndr.select(is_valid,mi.Float(1.),mi.Float(0.))deftraverse(self,callback:mi.TraversalCallback):# Registers the `g` parameter as a differentiable# parameter of the scenecallback.put('g',self._g,mi.ParamFlags.Differentiable)defto_string(self)->str:# Returns a humanly readable description of the materials=f"CustomRadioMaterial["\f"g={self._g}"\f"]"returns# We add a getter and setter to access `g`@propertydefg(self):returnself._g@g.setterdefg(self,v):self._g=mi.Float(v)# Register the custom radio material as a Mitsuba pluginmi.register_bsdf("custom-radio-material",lambdaprops:CustomRadioMaterial(props=props))
Let’s now use this custom radio material in a Sionna RT scene!We start by loading a new scene, then instantiate the newly created radiomaterial and use it for the reflector. We set\(g\) to a very low value toclearly see the impact on the reflected path gain.
scene_custom_mat=load_scene(rt.scene.simple_reflector,merge_shapes=False)scene_custom_mat.add(Transmitter("tx",position=(-2.,0.,1)))scene_custom_mat.add(Receiver("rx",position=(2.,0.,1)))scene_custom_mat.tx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")scene_custom_mat.rx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")my_mat=CustomRadioMaterial("custom-mat-instance",g=0.01)# Assign the custom radio material to the reflectorscene_custom_mat.objects["reflector"].radio_material=my_mat# To avoid confusion, discard the radio material initially loaded with the scenescene_custom_mat.remove("reflector-mat")# Print the material to visualize its parametersprint(scene_custom_mat.radio_materials)
{'custom-mat-instance':CustomRadioMaterial[g=[0.01]]}
We can see that theto_string() function wehave defined is used to describe the material.
We are now ready to trace paths.
paths=solver(scene_custom_mat)a,tau=paths.a,paths.taua=a[0].numpy()+1j*a[1].numpy()a=np.squeeze(a,(0,1,2,3))tau=tau.numpy()tau=np.squeeze(tau,(0,1))fora_,tau_inzip(a,tau):print("Delay [ns]: ",tau_*1e9," Gain [dB]",10.*np.log10(np.square(np.abs(a_))))
Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-76.33944
We can see that the gain of the reflected path is low because of the low value of\(g\).Let’s now set the value of\(g\) to a higher value:
my_mat.g=0.9# Paths coefficientpaths=solver(scene_custom_mat)a,tau=paths.a,paths.taua=a[0].numpy()+1j*a[1].numpy()a=np.squeeze(a,(0,1,2,3))tau=tau.numpy()tau=np.squeeze(tau,(0,1))fora_,tau_inzip(a,tau):print("Delay [ns]: ",tau_*1e9," Gain [dB]",10.*np.log10(np.square(np.abs(a_))))
Delay[ns]:13.342564Gain[dB]-55.370346Delay[ns]:14.917439Gain[dB]-56.79702
The reflected path now has a significantly higher gain!
Next, let’s compute a gradient of\(g\).The first step is to enable the gradient for this parameter.We then switch to theevaluated mode of as backpropagation throughadrjit loop can currently only be done with this mode.All we need to to next is to usedr.backward()to compute gradients.For this toy example, we use the total gain as objective function.
dr.enable_grad(my_mat.g)# Switch the computation of field loop to "evaluated" mode to# enable gradient backpropagation through the loopsolver.loop_mode="evaluated"paths=solver(scene_custom_mat)gain=dr.sum(dr.square(paths.a[0])+dr.square(paths.a[1]))dr.backward(gain)print("Gradient",my_mat.g.grad)
Gradient[2.32303e-06]
As expected, the gradient is positive as increasing the path gain requiresincreasing\(g\).
A More Complex Material Model
Let’s now enhance the previous radio material model by incorporating support forrefraction, which refers to radio waves passing through the material.Note that this model will not account for the refraction of radio waves as theypass through both the first and second boundaries of the scatterer.Instead, we will directly model the electromagnetic field that emerges afterpassing through the entire reflector.
In addition to the previous model for specular reflection described by(56)-(57),this enhanced model also supports refraction as follows.The transmitted wave\(\mathbf{E}_\text{t}\) propagates in the direction\(\mathbf{k}_\text{t}\) such that
where\(\hat{\mathbf{k}}_\text{i}\) is the direction of propagation of theincident wave.
The transmitted wave phasor\(\mathbf{E}_\text{t}\) is represented as
where
and
is the incident wave phasor.In these equations,\((\hat{\mathbf{e}}_{\text{i},s}, \hat{\mathbf{e}}_{\text{i},p})\)and\((\hat{\mathbf{e}}_{\text{t},s}, \hat{\mathbf{e}}_{\text{t},p})\) are the implicitbasis for the incident and scattered wave, respectively.
The next code snippet provides an implementation of the previous model as acustom radio material, with comments explaining every step.Note that, unlike the previous model, thesample()method is required to determine an interaction type (either specular reflection or refraction).This sampling process is executed by using a Bernoulli distribution where theprobability of selecting an interaction type corresponds to the ratio of energyscattered through that interaction.
Consulting the API documentation ofRadioMaterialBase is essentialfor a detailed understanding of this code.It is also recommended to consult theAPI documentation of the mi.BSDF class,asRadioMaterialBase inherits from it.
classEnhancedCustomRadioMaterial(RadioMaterialBase):# Note: The __init__ method is identical to the one of `CustomRadioMaterial`.# The __init__ method builds the radio material from:# - A unique `name` to identify the material instance in the scene# - The gain parameter `g`# - An optional `color` for displaying the material in the previewer and renderer# Providing these 3 parameters to __init__ is how an instance of this radio material# is built programmatically.## When loading a scene from an XML file, Mitsuba provides to __init__# only an `mi.Properties` object containing all the properties of the material# read from the XML scene file. Therefore, when a `props` object is provided,# the other parameters are ignored and should not be given.def__init__(self,name:str|None=None,g:float|mi.Float|None=None,color:Tuple[float,float,float]|None=None,props:mi.Properties|None=None):# If `props` is `None`, then one is built from the# other parametersifpropsisNone:props=mi.Properties("enhanced-custom-radio-material")# Name of the radio materialprops.set_id(name)props["g"]=gifcolorisnotNone:props["color"]=mi.ScalarColor3f(color)# Read the gain from `props`g=0.0if"g"inprops:g=props["g"]delprops["g"]self._g=mi.Float(g)# The other parameters (`name`, `color`) are given to the# base class to complete the initialization of the materialsuper().__init__(props)defsample(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,sample1:mi.Float,sample2:mi.Point2f,active:bool|mi.Bool=True)->Tuple[mi.BSDFSample3f,mi.Spectrum]:g=self._g# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Build the 3x3 change-of-basis matrix from the local basis to the world# basis.# `si.sh_frame` stores the three vectors that define the local interaction basis# in the world coordinate system.to_world=mi.Matrix3f(si.sh_frame.s,si.sh_frame.t,si.sh_frame.n).T# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# Direction of propagation of the transmitted field in the local coordinate systemkt_prime=ki_prime# Compute the Jones matrix in the implicit world coordinate system for the reflected# field and the transmitted field# The `jones_matrix_to_world_implicit()` builds the Jones matrix with the# structure we need.sqrt_g=mi.Complex2f(dr.sqrt(g),0.)jones_mat_ref=jones_matrix_to_world_implicit(c1=sqrt_g,c2=sqrt_g,to_world=to_world,k_in_local=ki_prime,k_out_local=kr_prime)sqrt_1mg=mi.Complex2f(dr.sqrt(1.-g),0.)jones_mat_tra=jones_matrix_to_world_implicit(c1=sqrt_1mg,c2=sqrt_1mg,to_world=to_world,k_in_local=ki_prime,k_out_local=kt_prime)# Sample the interaction type.# We use the `sample1` parameter, which is assumed to be a float uniformly sampled# from (0,1).# The probability of selecting a specular reflection corresponds to the ratio of energy# that is reflected, i.e., `g`.reflection=sample1<g# Select the Jones matrix and direction of the scattered wave according to the sampled# interaction typeko_prime=dr.select(reflection,kr_prime,kt_prime)jones_mat=dr.select(reflection,jones_mat_ref,jones_mat_tra)## We now only need to prepare the outputs# Cast the Jones matrix to a `mi.Spectrum` to meet the requirements of# the BSDF interface of Mitsubajones_mat=mi.Spectrum(jones_mat)# Instantiate and set the BSDFSample objectbs=mi.BSDFSample3f()# Specifies the type of interaction that was sampledbs.sampled_component=dr.select(reflection,InteractionType.SPECULAR,InteractionType.REFRACTION)# Direction of the scattered wave in the world framebs.wo=to_world@ko_prime# The next field of `bs` stores the probability that the sampled# interaction type and direction of scattering are sampled conditioned# on the given direction of incidence.bs.pdf=dr.select(reflection,g,1.-g)# Not used but required to be setbs.sampled_type=mi.UInt32(+mi.BSDFFlags.DeltaReflection)bs.eta=1.0returnbs,jones_matdefeval(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,wo:mi.Vector3f,active:bool|mi.Bool=True)->mi.Spectrum:g=self._g# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Build the 3x3 change-of-basis matrix from the local basis to the world# basis.# `si.sh_frame` stores the three vectors that define the local interaction basis# in the world coordinate system.to_world=mi.Matrix3f(si.sh_frame.s,si.sh_frame.t,si.sh_frame.n).T# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# Direction of propagation of the transmitted field in the local coordinate systemkt_prime=ki_prime# Read the sampled interaction type.# `si.prim_index` is used to store this information.sampled_event=si.prim_indexreflection=sampled_event==InteractionType.SPECULAR# Compute the Jones matrix in the implicit world coordinate system for the reflected# field and the transmitted field# The `jones_matrix_to_world_implicit()` builds the Jones matrix with the# structure we need.sqrt_g=mi.Complex2f(dr.sqrt(g),0.)jones_mat_ref=jones_matrix_to_world_implicit(c1=sqrt_g,c2=sqrt_g,to_world=to_world,k_in_local=ki_prime,k_out_local=kr_prime)sqrt_1mg=mi.Complex2f(dr.sqrt(1.-g),0.)jones_mat_tra=jones_matrix_to_world_implicit(c1=sqrt_1mg,c2=sqrt_1mg,to_world=to_world,k_in_local=ki_prime,k_out_local=kt_prime)# Select the Jones matrix and direction of the scattered wave according to the sampled# interaction typeko_prime=dr.select(reflection,kr_prime,kt_prime)jones_mat=dr.select(reflection,jones_mat_ref,jones_mat_tra)# This model only scatters energy in the direction of the specular reflection# and refraction.# Any other direction provided by the user `wo` should therefore lead to no energy.is_valid=isclose(dr.dot(ko_prime,wo),mi.Float(1.))jones_mat=dr.select(is_valid,jones_mat,0.)# Cast the Jones matrix to a `mi.Spectrum` to meet the requirements of# the BSDF interface of Mitsubajones_mat=mi.Spectrum(jones_mat)returnjones_matdefpdf(self,ctx:mi.BSDFContext,si:mi.SurfaceInteraction3f,wo:mi.Vector3f,active:bool|mi.Bool=True)->mi.Float:# Read the incident direction of propagation in the local coordinate# systemki_prime=si.wi# Direction of propagation of the reflected field in the local coordinate systemkr_prime=mi.reflect(-ki_prime)# Direction of propagation of the transmitted field in the local coordinate systemkt_prime=ki_prime# Read the sampled interaction type.# `si.prim_index` is used to store this information.sampled_event=si.prim_indexreflection=sampled_event==InteractionType.SPECULAR# Select the direction of the scattered wave and probability according to the sampled# interaction typeko_prime=dr.select(reflection,kr_prime,kt_prime)pdf=dr.select(reflection,self._g,1.-self._g)# The probability is set to `pdf` if the given direction `wo` matches# the expected direction for the scattered field `ko_prime`, or to# zero otherwiseis_valid=isclose(dr.dot(ko_prime,wo),mi.Float(1.))returndr.select(is_valid,pdf,mi.Float(0.))deftraverse(self,callback:mi.TraversalCallback):# Registers the `g` parameter as a differentiable# parameter of the scenecallback.put('g',self._g,mi.ParamFlags.Differentiable)defto_string(self)->str:# Returns a humanly readable description of the materials=f"EnhancedCustomRadioMaterial["\f"g={self._g}"\f"]"returns# We add a getter and setter to access `g`@propertydefg(self):returnself._g@g.setterdefg(self,v):self._g=mi.Float(v)# Register the custom radio material as a Mitsuba pluginmi.register_bsdf("enhanced-custom-radio-material",lambdaprops:EnhancedCustomRadioMaterial(props=props))
Note that we use thesample1 parameter ofsample()to select an interaction type.Thesample2 parameter is a 2-dimensional vector assumed to be uniformly sampledfrom the unit square\((0,1) \times (0,1)\)and used to sample a direction for the scattered wave when required (e.g.,when modeling diffuse reflections).
Let’s know use the enhanced custom material!We start by loading thesimple_reflector scene, but this time with weinstantiate two receivers: one on each side of the reflector to capture both areflected and a transmitted path.
scene_custom_mat=load_scene(rt.scene.simple_reflector,merge_shapes=False)scene_custom_mat.add(Transmitter("tx",position=(-2.,0.,1)))scene_custom_mat.add(Receiver("rx-1",position=(2.,0.5,1)))scene_custom_mat.add(Receiver("rx-2",position=(2.,-0.5,-1)))scene_custom_mat.tx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")scene_custom_mat.rx_array=PlanarArray(num_rows=1,num_cols=1,pattern="iso",polarization="V")# Instantiate newly created radio material and use it for the reflectormy_mat=EnhancedCustomRadioMaterial("custom-mat-instance",g=0.7)# Assign the custom radio material to the reflectorscene_custom_mat.objects["reflector"].radio_material=my_mat# To avoid confusion, discard the radio material initially loaded with the scenescene_custom_mat.remove("reflector-mat")# Print the material to visualize its parametersprint(scene_custom_mat.radio_materials)
{'custom-mat-instance':EnhancedCustomRadioMaterial[g=[0.7]]}
Tracing paths can be done as for the previous example.