Treatment of non-unit optic axis vectors

[gschlafly]

Description

When performing an iterative reconstruction, the optic axis estimates do not remain of unit length. This could be okay if only the direction of the vector was relevant in the forward model. However, the magnitude indirectly plays a role in the Jones calculus computations.

Files

https://github.com/PolarizedLightFieldMicroscopy/forward-model/blob/095e3eb406add7872218c7b549d77e4e6a2b72f5/VolumeRaytraceLFM/birefringence_implementations.py#L1102-L1106

To Reproduce

Functions extracted from birefringence_implementations.py:

import numpy as np
from VolumeRaytraceLFM.birefringence_implementations import  (
    JonesMatrixGenerators
)
def retardance(JM):
    '''Phase delay introduced between the fast and slow axis in a Jones Matrix'''
    e1, e2 = np.linalg.eigvals(JM)
    phase_diff = np.angle(e1) - np.angle(e2)
    retardance = np.abs(phase_diff)
    return retardance

def azimuth(JM):
    '''Rotation angle of the fast axis (neg phase)'''
    diag_sum = JM[0, 0] + JM[1, 1]
    diag_diff = JM[1, 1] - JM[0, 0]
    off_diag_sum = JM[0, 1] + JM[1, 0]
    a = np.imag(diag_diff / diag_sum)
    b = np.imag(off_diag_sum / diag_sum)
    azimuth = np.arctan2(-b, -a) / 2 + np.pi / 2
    return azimuth

def voxRayJM(Delta_n, opticAxis, rayDir, ell, wavelength):
    '''Compute Jones matrix associated with a particular ray and voxel combination'''
    # Azimuth is the angle of the slow axis of retardance.
    azim = np.arctan2(np.dot(opticAxis, rayDir[1]), np.dot(opticAxis, rayDir[2]))
    if Delta_n == 0:
        azim = 0
    elif Delta_n < 0:
        azim = azim + np.pi / 2
    ret = abs(Delta_n) * (1 - np.dot(opticAxis, rayDir[0]) ** 2) * 2 * np.pi * ell / wavelength
    JM = JonesMatrixGenerators.linear_retarder(ret, azim)
    return JM

The following example produces an incorrect retardance and azimuth angle.

opticAxis = np.array([1,1,0])
JM = voxRayJM(0.01, opticAxis, [np.array([1,0,0]), np.array([0,1,0]), np.array([0,0,1])], 1, 0.5)
ret = retardance(JM)
azim = azimuth(JM)
print(f"ret: {ret}, azim: {azim}")

Solution ideas

1. Optic axis can be normalized before performing the Jones calculus operations. For pytorch gradient purposes, we would need to define a new variable.
2. When taking the dot product, we can divide by the norm of the optic axis.
3. We can parameterize the optic axis vector as two angles instead of a 3D vector.

[gschlafly]

51 includes investigations about the norm of the optic axis.

[gschlafly]

With the assert statements added in commit 2e75ba82b8944b95087f2f26c11c54e95c27b7ea to the estimated-vol branch to check that the optic axis vector do not have norm zero, 9 of the 42 tests fail.

[gschlafly]

Similar issue #61:

If the optic axis is not a unit vector, the projection of the ray onto the optic axis is misleading. In e4a27d7, we normalized the optic axis in the projection. The same should be done for the torch process.

[gschlafly]

Optic axis is normalized after each iteration update.