5.2.4. Analysis of experimental Mueller Matrix

In the first example, we saw how we can measure the Mueller matrix of an unknowns sample. This one corresponds to the second part: extract information from that Mueller matrix. This is implemented in the analysis, checks and parameters methods.

[1]:
from py_pol import np, degrees
from py_pol.mueller import Mueller
import matplotlib.pyplot as plt
3.0.2

5.2.4.1. Filtering

Experimental measurements are subject to errors that distorts the Mueller matrix recovered from a real polarimeter. Fortunately, there are some mathematical tools we can use in order to reduce them. This is usually refered as filtering.

5.2.4.1.1. Physical conditions

This type of filtering is a must for all experimental measurements. It is based on the physical properties that Mueller matrices of real objects must satisfy. With the method filter_physical_conditions, we transform our Mueller matrix in such a way that it fulfills all those conditions. The user may consult the documentation of the method is_physical (Check_Mueller class) for more information about those conditions.

Example: Experimental matrix of air.

[2]:
# Create The object of the experimental and theorical matrices
m_air = np.matrix([[0.9924, -0.0238, -0.0096, -0.0024],
              [0.0198, 1.0047, -0.0066, 0.0004],
               [0.0058, 0.0055, 1.0086, 0.0050],
               [-0.0003, 0.0053, 0.0011, 0.9924]
              ]) / 0.9924
M_exp = Mueller('Experimental')
M_exp.from_matrix(m_air)
M_target = Mueller('Theor')
M_target.vacuum()

# Filter
M_filtered = M_exp.analysis.filter_physical_conditions(tol=1e-2, verbose=True, keep=True)
M_filtered.name = 'Filtered'
The original matrix is:
[[ 1.00000000e+00 -2.39822652e-02 -9.67351874e-03 -2.41837969e-03]
 [ 1.99516324e-02  1.01239420e+00 -6.65054414e-03  4.03063281e-04]
 [ 5.84441757e-03  5.54212011e-03  1.01632406e+00  5.03829101e-03]
 [-3.02297461e-04  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
Second condition was violated. Fixed matrix is:
[[ 1.00000000e+00 -2.39822652e-02 -9.67351874e-03 -2.41837969e-03]
 [ 1.99516324e-02  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 5.84441757e-03  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [-3.02297461e-04  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
The original matrix is:
[[ 1.00000000e+00 -2.39822652e-02 -9.67351874e-03 -2.41837969e-03]
 [ 1.99516324e-02  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 5.84441757e-03  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [-3.02297461e-04  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
Fifth A condition was violated. Fixed matrix is:
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 1.99516324e-02  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 5.84441757e-03  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [-3.02297461e-04  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
The original matrix is:
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 1.99516324e-02  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 5.84441757e-03  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [-3.02297461e-04  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
Fifth B condition was violated. Fixed matrix is:
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 0.00000000e+00  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [ 0.00000000e+00  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
The original matrix is:
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  1.00000000e+00 -6.65054414e-03  4.03063281e-04]
 [ 0.00000000e+00  5.54212011e-03  1.00000000e+00  5.03829101e-03]
 [ 0.00000000e+00  5.34058847e-03  1.10842402e-03  1.00000000e+00]]
None condition was violated.

We can check if the filtering procedure has improved or not the experimental Mueller matrix.

[3]:
# Calculate mean square error from the ideal case
error_exp = np.linalg.norm(M_target.M - M_exp.M)/16
error_filter = np.linalg.norm(M_target.M - M_filtered.M)/16
# Print results
print("The experimental matrix is:")
print(M_exp)
print("- MSE of the experimental matrix: {}.\n\n".format(error_exp))
print("The filtered matrix is:")
print(M_filtered)
print("- MSE of the filtered matrix: {}.".format(error_filter))
The experimental matrix is:
Experimental =
  [+1.0000, -0.0240, -0.0097, -0.0024]
  [+0.0200, +1.0124, -0.0067, +0.0004]
  [+0.0058, +0.0055, +1.0163, +0.0050]
  [-0.0003, +0.0053, +0.0011, +1.0000]

- MSE of the experimental matrix: 0.00254431357763787.


The filtered matrix is:
Filtered =
  [+1.0000, +0.0000, +0.0000, +0.0000]
  [+0.0000, +1.0000, -0.0067, +0.0004]
  [+0.0000, +0.0055, +1.0000, +0.0050]
  [+0.0000, +0.0053, +0.0011, +1.0000]

- MSE of the filtered matrix: 0.0007132730401082221.

The error has been reduced in a factor of ~4. We can also see that the elemnts of the diagonal are almost exactly 1, and the errors in the rest of the ellementrs are at most in the 3rd decimal.

5.2.4.1.2. Purify

Many times, we know beforehand some information about the optical element that is measured by a polarimeter. One of the most common facts we can know is if the element does not depolarize, i.e., it is a pure element. Due to the experimental errors, most probably the Mueller matrix retrieved by the polarimeter will not be pure.

A pure Mueller matrix is pure if all the eigenvalues of the covariance matrix are 0 except one. We have implemented two methods that works similarly in order to purify the Mueller matrix: filter_purify_number and filter_purify_threshold. Both make 0 some eigenvalues of the covariance matrix.

Example: Transmission measurement of a chiral turbid famtom (Ghosh, N. et al. 2008. Mueller matrix decomposition for extraction of individual polarization parameters from complex turbid media exhibiting multiple scattering, optical activity, and linear birefringence. Journal of Biomedical Optics 13, 044036-14.). This example is performed mathematically in J.J. Gil, R. Ossikovsky “Polarized light and the Mueller Matrix approach”, CRC Press (2016) pp 226.

[4]:
from py_pol.utils import order_eig

# Create the object of the experimental measurement
m_exp = np.matrix([[1, -0.0229, 0.0027, 0.0058],
              [-0.0186, 0.9956, -0.0361, 0.0318],
               [-0.0129, 0.0392, 0.2207, -0.9656],
               [0.0014, 0.028, 0.9706, 0.2231]
              ]) # This is a normalized matrix, so we will always show normalized ones
M_exp = Mueller('Experimental')
M_exp.from_matrix(m_exp)

# Calculate the covariance matrix
H = M_exp.covariance_matrix()
print("The covariance matrix is:")
print(H)

# Calculate eigenvalues
qi, U = np.linalg.eig(H)
qi, U = order_eig(qi, U)
print("\nThe eigenvalues are:")
print(qi)
The covariance matrix is:
[[ 4.88525e-01+0.j      -8.35000e-03+0.0094j   6.57500e-03-0.00735j
   1.10950e-01-0.48405j]
 [-8.35000e-03-0.0094j   2.17500e-03+0.j      -6.00000e-04-0.00125j
  -1.30250e-02+0.00665j]
 [ 6.57500e-03+0.00735j -6.00000e-04+0.00125j  2.50000e-05+0.j
   9.70000e-03-0.0065j ]
 [ 1.10950e-01+0.48405j -1.30250e-02-0.00665j  9.70000e-03+0.0065j
   5.09275e-01+0.j     ]]

The eigenvalues are:
[ 9.96216864e-01  4.02577388e-03  8.03161335e-04 -1.04579885e-03]

The last eigenvalue is lower than 0, which is not true in physical Mueller matrices. It must be made zero so M is physical (this operation is also performed by filter_physical_conditions, along with other ones).

[5]:
M_filter1 = M_exp.analysis.filter_purify_number(Neig=1, keep=True)
M_filter1.name = 'Filter 1 eig'
print(M_filter1/M_filter1.m00)
Filter 1 eig =
  [+1.0000, -0.0222, +0.0023, +0.0059]
  [-0.0193, +0.9937, -0.0358, +0.0314]
  [-0.0127, +0.0396, +0.2204, -0.9651]
  [+0.0017, +0.0281, +0.9689, +0.2229]

d:\codigo\py_pol\py_pol\mueller.py:532: ComplexWarning: Casting complex values to real discards the imaginary part
  M[i, j] = elem

As the second and third eigenvalues are very small, it is reasonable to think that the measured optical element is pure. So we can filter all three eigenvalues.

[6]:
M_filter3 = M_exp.analysis.filter_purify_number(Neig=3, keep=True)
M_filter3.name = 'Filter 3 eig'
print(M_filter3/M_filter3.m00)

cond = M_filter3.checks.is_pure()
print("The matrix M_filter3 is pure? Answer: {}".format(cond))
Filter 3 eig =
  [+1.0000, -0.0211, +0.0009, +0.0091]
  [-0.0208, +0.9988, -0.0362, +0.0320]
  [-0.0095, +0.0395, +0.2226, -0.9739]
  [+0.0024, +0.0281, +0.9740, +0.2238]

The matrix M_filter3 is pure? Answer: True

5.2.4.2. Behavior analysis

The first thing that you may do in order to try to understand the meaning of a Mueller matrix is analyze which is its behavior as diattenuator or retarder. There are three methods that can be used: diattenuator, polarizer and retarder. Those methods analyze the mueller matrix AS IF they correspond to an homogeneous diattenuator, a polarizer or a retarder. Take a note that it is a big if.

[7]:
M1 = Mueller('diattenuator')
M1.diattenuator_charac_angles_from_vector(p1=0.75, p2=0.25, alpha=25*degrees, delay=190*degrees)
print(M1)

p1, p2, alpha, delay = M1.analysis.diattenuator(param='charac')
print("M1 has the following diattenuator parameters:")
print("  Transmission coefs: -p1 = {}.           -p2 = {}.".format(p1, p2))
print("  Transparent state:  -alpha = {} deg.    -delay = {} deg.".format(alpha/degrees, delay/degrees))
diattenuator =
  [+0.3125, +0.1607, -0.1886, -0.0333]
  [+0.1607, +0.2391, -0.0606, -0.0107]
  [-0.1886, -0.0606, +0.2586, +0.0125]
  [-0.0333, -0.0107, +0.0125, +0.1897]

M1 has the following diattenuator parameters:
  Transmission coefs: -p1 = 0.75.           -p2 = 0.24999999999999997.
  Transparent state:  -alpha = 25.0 deg.    -delay = 190.0 deg.

As the Mueller matrix has no errors, the measurement of parameters is exact.

But let’s see what happens in a case where we have errors. We can use these methods to check the usefullness of the filtering we saw in the previous section. From the nature of the sample, we can know that the sample must behave approximately as a retarder.

[8]:
# Experimental matrix
_, errors = M_exp.checks.is_retarder(give_all=True)
retardance, azimuth, ellipticity = M_exp.analysis.retarder(param='azel')
print("Distance from Experimental matrix to be considered a pure retarder: {}.".format(
    np.linalg.norm(errors)))
print("                   - Retardance = {} deg.".format(retardance/degrees))
print("  Fast eigenstate: - Azimuth           = {} deg".format(azimuth/degrees))
print("                   - Ellipticity angle = {} deg".format(ellipticity/degrees))
Distance from Experimental matrix to be considered a pure retarder: 0.02289734994415676.
                   - Retardance = 77.3085868316892 deg.
  Fast eigenstate: - Azimuth           = 90.0562244816153 deg
                   - Ellipticity angle = -1.1135705865327494 deg
[9]:
# Filtering of one eigenvalue
_, errors = M_filter1.checks.is_retarder(give_all=True)
retardance, azimuth, ellipticity = M_filter1.analysis.retarder(param='azel')
print("Distance from Filter 1 matrix to be considered a pure retarder: {}.".format(
    np.linalg.norm(errors)))
print("                   - Retardance = {} deg.".format(retardance/degrees))
print("  Fast eigenstate: - Azimuth           = {} deg".format(azimuth/degrees))
print("                   - Ellipticity angle = {} deg".format(ellipticity/degrees))
Distance from Filter 1 matrix to be considered a pure retarder: 0.0262939870146039.
                   - Retardance = 77.37974086749291 deg.
  Fast eigenstate: - Azimuth           = 90.04951588351138 deg
                   - Ellipticity angle = -1.1163055731148146 deg
[10]:
# Filtering of 3 eigenvalues
_, errors = M_filter3.checks.is_retarder(give_all=True)
retardance, azimuth, ellipticity = M_filter3.analysis.retarder(param='azel')
print("Distance from Filter 1 matrix to be considered a pure retarder: {}.".format(
    np.linalg.norm(errors)))
print("                   - Retardance = {} deg.".format(retardance/degrees))
print("  Fast eigenstate: - Azimuth           = {} deg".format(azimuth/degrees))
print("                   - Ellipticity angle = {} deg".format(ellipticity/degrees))
Distance from Filter 1 matrix to be considered a pure retarder: 0.010861365008514036.
                   - Retardance = 77.13732793253152 deg.
  Fast eigenstate: - Azimuth           = 90.05750862281892 deg
                   - Ellipticity angle = -1.1134522341557618 deg

In this case, the filtering does not affect much in the measured parameters. But we can see that after filtering, we are much more certain that the element behaves as a pure retarder.

5.2.4.3. Decomposition

Using the methods described in the previous section, we need to know beforehand if the element we want to analyze corresponds to a ddiattenuator, polarizer or retarder. That is not always the case. Also, it is usual that real optical elements show diattenuation, polarizance, retardance and depolarization at the same time.

Fortunately, due to their properties, all physical Mueller matrices can be decomposed into a depolarizer, an homogenous retarder and an homogenous diattenuator. We can use this polar decomposition in order to understand the properties of an optical element.

Example: chiral turbid famtom from previous sections

[11]:
Md, Mr, Mp = M_exp.analysis.decompose_polar(verbose=True)
------------------------------------------------------
General case.
Polar decomposition of the matrix M = Mdesp * Mr * Mp:

The depolarizer Mueller matrix is:
depolarizer =
  [+1.0000, +0.0000, +0.0000, +0.0000]
  [+0.0041, +0.9969, +0.0000, -0.0000]
  [-0.0070, +0.0000, +0.9915, -0.0000]
  [-0.0019, -0.0000, -0.0000, +0.9966]

Parameters:
  - Polarizance = 0.008334694965594261.
  - Depolarization degree = 0.09985338379915462.

The retarder Mueller matrix is:
retarder =
  [+1.0000, +0.0000, +0.0000, +0.0000]
  [+0.0000, +0.9988, -0.0362, +0.0320]
  [+0.0000, +0.0393, +0.2227, -0.9741]
  [+0.0000, +0.0281, +0.9742, +0.2239]

Parameters:
  - Delay = 77.1309632136399 deg.
  - Angle = 88.8894601282744 deg; Delay between components = 267.0616037216929 deg.
  - Azimuth = 90.05695718867285 deg; Ellipticity = -1.109079035173538 deg.

The diatenuator/polarizer Mueller matrix is:
diattenuator =
  [+1.0000, -0.0229, +0.0027, +0.0058]
  [-0.0229, +1.0000, -0.0000, -0.0001]
  [+0.0027, -0.0000, +0.9997, +0.0000]
  [+0.0058, -0.0001, +0.0000, +0.9997]

Parameters:
  - p1 = 1.011818600108798; p2 = 0.9880400399143103.
  - Angle = 82.19551946316135 deg; Delay between components = 65.03721016919296 deg.
  - Azimuth = 86.63782345344518 deg; Ellipticity = 7.059435856767925 deg.

The mean square error in the decomposition is: 1.1789986919039411e-15
The maximum error in the decomposition is: 9.992007221626409e-16
------------------------------------------------------

From this decomposition, we can see that the sample barely depolarizes (Depolarization degree of the depolarizer = 0.10). Also, the difference between p1 and p2 is small (1.01 and 0.99 respectively). The parameters of the retarder are very similar to the ones determined in the previous section, but their accuracy will be higher, as the depolarization, polarizance and diattenuation have been isolated from retardance.

Note: p1 is higher than 1. This is due to the fact that M_exp is not a physical Mueller matrix.

Now, we can do the same with the filtered matrix. As the matrix has been forced to be pure, we can use the method decompose_pure instead of decompose_polar, which decomposes any pure Mueller matrix into a diattenuator and a retarder.

[16]:
Mr, Md = M_filter3.analysis.decompose_pure(verbose=True)
------------------------------------------------------
Matrx M decomposed as M = Mr * Md.

The retarder Mueller matrix is:
retarder =
  [+1.0000, -0.0000, +0.0000, +0.0000]
  [+0.0000, +0.9988, -0.0362, +0.0321]
  [-0.0000, +0.0394, +0.2227, -0.9741]
  [+0.0000, +0.0281, +0.9742, +0.2239]

Parameters:
  - Delay = 77.13274318504531 deg.
  - Angle = 88.88659355522645 deg; Delay between components = 266.98801477841516 deg.
  - Azimuth = 90.05853316773448 deg; Ellipticity = -1.1118675730541951 deg.

The diatenuator Mueller matrix is:
diattenuator =
  [+0.9962, -0.0210, +0.0009, +0.0090]
  [-0.0210, +0.9962, -0.0000, -0.0001]
  [+0.0009, -0.0000, +0.9960, +0.0000]
  [+0.0090, -0.0001, +0.0000, +0.9960]

Parameters:
  - p1 = 1.0095126177435985; p2 = 0.9865688024040568.
  - Angle = 78.31151140019138 deg; Delay between components = 84.06091902472883 deg.
  - Azimuth = 88.71953234482453 deg; Ellipticity = 11.622053674000005 deg.

The mean square error in the decomposition is: 1.2762915567195726e-16
The maximum error in the decomposition is: 1.1102230246251565e-16
------------------------------------------------------