As Alan has already said, following the actual path of the light raypath through each layer leads to more physically accurate results. I will base my answer on a paper by Andrea Weidlich and Alexander Wilkie ("Arbitrarily Layered Micro-Facet Surfaces") that I have read and partially implemented. In theirthe paper theythe authors assume that the distance between two layers is smaller than the radius of a differential area element. This assumption simplifies the implementation because we do not have to calculate intersection points separately for each layer, actually we assume that theall intersection points over the layers are just the same for all layerspoint.
UPDATE: Actually I have adopted a different method to implement the evaluation of this layered model. While I have kept with the idea of considering the intersection points to be just the same point along the layers, I have computed the sampling and the final BRDF differently: for sampling, I have used ordinary ray tracing, but through the layers (using Russian Roulette to select between reflection/refraction when that's the case); for the final BRDF evaluation, I just multiply each BRDF traversed by the ray path (weighting the incident radiances according to the cosine of the incident ray).
In this first stage we will determine the actual light path through the layers. When a light ray is moving from a less dense medium, e.g. air, to a more densedenser medium, e.g. glass, part of its energy is reflected and the remaining part is transmitted. You can find the amount of energy that is reflected through the Fresnel reflectance equations. So, for instance, if the Fresnel reflectance of a given dielectric is 0.3, we know that 30% of the energy is reflected and 70% will be transmitted:
When the light ray is moving from a more densedenser to a less dense medium, the same principle described by the Fresnel reflectance applies. However, in this specific case, total internal reflection (a.k.a TIR) might also happen if the angle of the incident ray is above the critical angle. In the case of TIR, 100% of the energy is reflected back into the material:
When light hits a conductor or a diffuse surface, it will always be reflected (being the direction of reflection related to the type of the BRDF). In a multilayer material, the resulting light path will be the agregateaggregate result of all those possibilities. Thus, in the case of a 3-layer material, assuming that the first and secongsecond layers are dielectrics and the third layer is diffuse, we might end up, for instance, with the following light path (a tree actually):
We can simulate this type of interaction using recursion and weighting each light path according to the actual reflectance/transmitance at the corresponding incident pointpoints. A problem regarding the use of recursion in this case is that the number of rays increases with the deepness of the recursion, concentrating computational effort on rays that individually might contribute almost nothing to the final result. On the other hand, the aggregate result of those individual rays at deep recursion levels can be significant and should not be discarded. In this case, we can use Russian Roulette (RR) in order to avoid branching and to probabilistic end light paths without losing energy, but at the cost of a higher variance (noisier result). In this case, the result of the Fresnel reflectance, or the TIR, will be used to randomly select which path to follow. For instance:
As can be seen, TIR or Fresnel reflectance might keep some rays bouncing indefinitely among layers. As far as I know, Mitsuba implements plastic as a two layer material, and it uses a closed form solution for this specific case that accounts for an infinity number of light bounces among layers. However, Mitsuba also allows for the creation of multilayer materials with an arbitrary number of layers, in which case it imposes a maximum number of internal bounces since no closed form solution seems to exist for the general case. As a side effect, some energy can be lost in the rendering process, making the material look darker than it should eventually be.
An additional option is to mix branching and RR. For instance, the initial rays (lower deep levels) might present substantial contribution to the final image. Thus, one might choose to branch only at the first one or two intersections, using only RR afterwards. This is, for example, the case withapproached used by smallpt.
An interesting point regarding multilayered materials is that individual reflected/transmitted rays can be importance sampled according to the corresponding BRDFs/BTDFs of the currenteach layer.
Since all terms of the estimator are multiplied, we can simplify the implementation by computing the final BSDF and $pdf$ and factoring out the $L_i$ term:
$$fr = fr_1 \cdot fr_2 \cdot fr_3 \cdot fr_2 \cdot fr_1$$
$$pdf = pdf_1 \cdot pdf_2 \cdot pdf_3 \cdot pdf_2 \cdot pdf_1$$
$$\cos \theta= \cos \theta_1 \cdot \cos \theta_2 \cdot \cos \theta_3 \cdot \cos \theta_2 \cdot \cos \theta_1$$
$$ L_r = \left( \frac{fr \cos \theta}{pdf} \right) L_i$$
The paper by Andrea Weidlich and Alexander Wilkie also takes absorption into consideration, i.e. each light ray might be attenuated according to the absorption factor of each transmissive layer and to the distance traveled by the ray within the layerits thickness. I've not included absorption into my renderer yet, but it is represented by just a real coefficient computedone scalar value, which will be evaluated according to the Beer's Law.