Derivation of delta_x and delta_y error bounds for triangle Intersection

[Erfan-Ahmadi]

https://github.com/mmp/pbrt-v3/blob/aaa552a4b9cbf9dccb71450f47b268e0ed6370e2/src/shapes/triangle.cpp#L276-L280

As described in the book 3.9 and 3.6 there are three transformations applied to the triangle to continue with intersection test in the new coordinate system.

Calculation of a new_x value for a vertex looks like this equation:

new_x = (p0_x - Origin_x) + (-dx * (1/dz)) * (p0_z - Origin_z)

all the operations in the above are floating point and after error analysis and initial simplification to gamma functions we arrive at:

new_x ⊂ (1+-gamma_2) * (p0_x - Origin_x) + (1+-gamma_5) * (-dx) * (p0_z - Origin_z) * (1/dz)

I cannot derive the final formula after simplifications of the bounds in formula above.

It is also surprising to me that no explanation is given in the book or github code.

I would like to know if the final formula is derived from simplifing the above equation and what assumptions were used to simplify it to:

delta_x = gamma_5 * ( Max | Xi | + Max | Zi | )

I know that:

| (p0_z - Origin_z) * (1/dz) | <= Max |Zi| and |(p0_x-Origin_x) - (dx / dz) * (p0_z - Origin_z)| <= Max |Xi|

I have derived all the error bounds formulas in the book myself but this one I cannot figure out.

(note that Xi and Zi are the transformed vertices components and are NOT the initial X and Z's for the mesh)

[kennyalive]

I checked the error bounds derivation for x,y coordinates in "Watertight Ray/Triangle Intersection" article by Sven Woop, et al., section 3.3 and it looks good for me. They used (1+x)^n -> (1 + (n+1)*x) approximation and they also got machine epsilon multiplied by 5 but different other constants. Also if you use gamma() approximation in their approach it will result in tighter bounds (machine epsilon * 4)

UPDATE: It looks I got similar error bounds as in the book with a difference that it uses translated coordinates instead of final translation+shear coordinates and gamma_4 instead of gamma_5:

delta_x = gamma_4 * (Max(Xi) + Max(Zi)) Xi = Xi_original - ray.origin Zi = Zi_original - ray.origin

At first I run error analysis strictly according to computations of the algorithm, then I had epression similar to gamma2 * a + gamma4 * S * b where S is in [0..1]. This expression can be simplifying by replacing S with its max value - one, and also by replacing gamma2 with gamma4. This gets less tight but still conservative bounds+simple expression: gamma4*(a + b)