{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using the Mueller module \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we demonstrate how to use the Mueller module, which defines the `MuellerMatrix` and `StokesVector` classes and various operations on them.\n", "\n", "The first thing, of course, is to import the package and some other packages that we will be using. We will be creating some random matrices and vectors, so we need `random`, and we will be graphing a few things, so we need `maptplotlib.pyplot`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pySCATMECH.mueller import *\n", "import random \n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mueller matrices\n", "\n", "The class `MuellerMatrix` is used for handling Mueller matrices. Here we show several ways one can create a unit Mueller matrix. The `JonesMueller()` function converts a 2 x 2 complex Jones matrix into a Mueller matrix. And `MuellerUnit()` returns a unit Mueller matrix." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0. 0. 0.]\n", " [0. 1. 0. 0.]\n", " [0. 0. 1. 0.]\n", " [0. 0. 0. 1.]]\n" ] } ], "source": [ "m1 = MuellerMatrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])\n", "m2 = MuellerMatrix([1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1])\n", "m3 = MuellerMatrix(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1)\n", "m4 = JonesMueller([[1,0],[0,1]])\n", "m5 = MuellerUnit()\n", "print(m5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check to see that these are valid matrices. The `valid` method tests whether a Mueller matrix maps the set of valid Stokes vectors onto itself. The `physically_valid` method tests whether a Mueller matrix is a convex sum of Jones-Mueller matrices. The latter is a more stringent test of whether the matrix is realizable." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "True\n" ] } ], "source": [ "print(m1.valid())\n", "print(m1.physically_valid())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But, a random matrix is not necessarily valid. In fact, it probably isn't!" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number valid() is True: 1\n", "Number physically valid() is True: 0\n" ] } ], "source": [ "# Returns a random number with a normal distribution...\n", "ran = lambda : random.gauss(0,1)\n", "\n", "# Create 20 random matrices...\n", "m1 = [MuellerMatrix([[ran(), ran(), ran(), ran()],\n", " [ran(), ran(), ran(), ran()],\n", " [ran(), ran(), ran(), ran()],\n", " [ran(), ran(), ran(), ran()]]) for i in range(20)]\n", "\n", "# See if they are valid...\n", "print(\"Number valid() is True:\", sum([m.valid() for m in m1]))\n", "print(\"Number physically valid() is True:\", sum([m.physically_valid() for m in m1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A Mueller matrix made from a random Jones matrix, however, always yields a valid Mueller matrix." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number value() is True: 20\n", "Number physically valid() is True: 20\n" ] } ], "source": [ "# Returns a random complex number with a normal distribution...\n", "ranc = lambda : random.gauss(0,1)+random.gauss(0,1)*1j\n", "\n", "# Create 20 random Jones-Mueller matrices...\n", "m1 = [JonesMueller([[ranc(), ranc()], [ranc(), ranc()]]) for i in range(20)]\n", "\n", "# See if they are valid...\n", "print(\"Number value() is True:\", sum([m.valid() for m in m1]))\n", "print(\"Number physically valid() is True:\", sum([m.physically_valid() for m in m1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make a random depolarizing Mueller matrix, we average some random Jones-Mueller matrices:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0.79126519 0.51754966 0.04396498 -0.09197144]\n", " [-0.21857972 -0.10891561 0.03399935 0.04409256]\n", " [ 0.24109773 0.27199561 0.06684414 -0.02241909]\n", " [ 0.31260027 0.21946993 0.16903187 -0.02210762]]\n", "m is a valid Stokes to Stokes mapper: True\n", "m is the sum of Jones-Mueller matrices: True\n" ] } ], "source": [ "def RandomMuellerMatrix(n):\n", " \"\"\"\n", " Creates a random Mueller matrix by averaging `n` Jones-Mueller matrices.\n", " \"\"\"\n", " ranc = lambda : (random.gauss(0,1)+random.gauss(0,1)*1j)/2\n", " return sum([JonesMueller([[ranc(),ranc()],[ranc(),ranc()]]) for i in range(n)])/n\n", "\n", "m = RandomMuellerMatrix(4)\n", "print(m)\n", "print(\"m is a valid Stokes to Stokes mapper:\",m.valid())\n", "print(\"m is the sum of Jones-Mueller matrices:\",m.physically_valid())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stokes Vectors\n", "\n", "\n", "The class `StokesVector` is used for handling Stokes vectors.\n", "We can create a Stokes vector a several ways. The `Polarization` function provides a means for generating arbitrary polarizations by name or parameter. (Note that SCATMECH's convention for the signs of the Stokes vector elements differs from some others.)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Unpolarized Stokes vector: StokesVector(1,0,0,0)\n", "Unpolarized Stokes vector: StokesVector(1,0,0,0)\n", "s-polarized unit polarization from a Jones vector: StokesVector(1.0,1.0,0.0,0.0)\n", "p-polarized unit Stokes vector: StokesVector(1,-1,0,0)\n", "s-polarized unit Stokes vector: StokesVector(1,1,0,0)\n", "Left-handed circularly polarized unit Stokes vector: StokesVector(1,0,0,1)\n", "Right-handed circularly polarized unit Stokes vector with intensity 2: StokesVector(2,0,0,-2)\n", "Unpolarized Stokes vector: StokesVector(1,0,0,0)\n", "Unpolarized Stokes vector: StokesVector(1,0,0,0)\n", "Linear polarization at 45 degrees: StokesVector(1.0,6.123233995736766e-17,1.0,0.0)\n", "Elliptically polarized radiation: StokesVector(1.0,0.2869652992419897,0.7102640395405222,0.6427876096865393)\n", "Partially polarized radiation: StokesVector(1.0,0.8,0.0,0.0)\n", "Partially polarized elliptical polarized radiation with intensity 3:\n", " StokesVector(3.0,0.6026271284081782,1.4915544830350964,1.3498539803417322)\n" ] } ], "source": [ "print(\"Unpolarized Stokes vector:\", StokesVector((1,0,0,0)) )\n", "\n", "print(\"Unpolarized Stokes vector:\", StokesVector(1,0,0,0) )\n", "\n", "print(\"s-polarized unit polarization from a Jones vector:\", JonesStokes((1.,0.)) )\n", "\n", "print(\"p-polarized unit Stokes vector:\", Polarization('p') )\n", "\n", "print(\"s-polarized unit Stokes vector:\", Polarization(state = 's') )\n", "\n", "print(\"Left-handed circularly polarized unit Stokes vector:\", Polarization(\"L\") )\n", "\n", "print(\"Right-handed circularly polarized unit Stokes vector with intensity 2:\", Polarization(\"R\", I=2) )\n", "\n", "print(\"Unpolarized Stokes vector:\", Polarization(\"U\") )\n", "\n", "print(\"Unpolarized Stokes vector:\", Polarization('s', DOP=0) )\n", "\n", "print(\"Linear polarization at 45 degrees:\", Polarization(angle = 45*deg) )\n", "\n", "print(\"Elliptically polarized radiation:\", Polarization(angle = 34*deg, ellipticity = 20*deg) )\n", "\n", "print(\"Partially polarized radiation:\", Polarization(\"S\", DOP = 0.8) )\n", "\n", "print(\"Partially polarized elliptical polarized radiation with intensity 3:\\n\", \n", " Polarization(angle = 34*deg, ellipticity = 20*deg, DOP = 0.7, I = 3) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that polarization sensitivity can also be expressed as a Stokes vector. However, there is often a factor of two that enters for all but unpolarized sensitivity. That is, a perfect polarizer transmits 50 % of unpolarized radiation. We can apply `senstivity = True` to the `Polarization()` function, or use a `Sensitivity` function." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Signal for p-polarization: 1.0\n", "Signal for s-polarization: 0.0\n" ] } ], "source": [ "# Only sensitive to p-polarization\n", "# The following two statements are equivalent...\n", "sens = Polarization('p', sensitivity = True)\n", "sens = Sensitivity('p')\n", "\n", "print('Signal for p-polarization:',sens @ Polarization('p'))\n", "print('Signal for s-polarization:',sens @ Polarization('s'))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sensitivity does not need to be full polarizing. By setting `DOP` in the range [0,1], sensitivity to the orthogonal polarization can be non-zero." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sens = StokesVector(0.6,-0.4,0.0,0.0)\n", "Signal for p-polarization: 1.0\n", "Signal for s-polarization: 0.19999999999999996\n" ] } ], "source": [ "sens = Sensitivity('p',DOP=.8)\n", "print('sens = ', sens)\n", "\n", "print('Signal for p-polarization:',sens @ Polarization('p'))\n", "print('Signal for s-polarization:',sens @ Polarization('s'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a manner similar to what we did with Mueller matrices, we can create a random Stokes vector by averaging a number of Stokes \n", "vectors created from Jones vectors. For a single Jones vector, the Stokes vector is depolarized. For more than one, the Stokes vector is partially polarized or unpolarized. We illustrate that by plotting the degree of polarization (using the `DOP()` method) as a function of the number of Jones-Stokes vectors we add." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-03-02T13:10:59.362924\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.5.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def RandomStokesVector(n):\n", " ranc = lambda : (random.gauss(0,1) + random.gauss(0,1)*1j)/2\n", " return sum([JonesStokes([ranc(),ranc()]) for i in range(n)])/n\n", "\n", "plt.figure()\n", "plt.plot(range(1,100),[RandomStokesVector(i).DOP() for i in range(1,100)])\n", "plt.title(\"Degree of polarization vs. number of Jones-Stokes vectors\")\n", "plt.xlabel(\"n\")\n", "plt.ylabel(\"Degree of polarization\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can multiply Mueller matrices and Stokes vectors. You can use the `@` or `*` operators or the `dot()` function to multiply the objects. If a Stokes vector multiplies a Mueller matrix from the left, it is assumed to be transposed first. The product of two Stokes vectors is the inner product, a scalar. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "m1 = \n", " [[ 0.94335668 -0.16250793 0.46915549 -0.14212114]\n", " [-0.23542067 0.62314123 -0.18648852 -0.05825565]\n", " [-0.01326436 -0.06702373 0.02740091 -0.4182326 ]\n", " [ 0.37937625 -0.06489837 0.54109253 0.07621696]]\n", "m2 = \n", " [[ 0.90393586 0.42883859 0.05893816 -0.0398583 ]\n", " [-0.03255312 -0.05484299 0.28021642 0.20824843]\n", " [ 0.17437584 0.3752689 -0.03803409 0.08744445]\n", " [ 0.35288836 0.34723155 0.26526746 -0.43082467]]\n", "s1 = StokesVector(1.2205852882966242,-0.07508395251618234,0.8002605801487375,-0.2529729329906206)\n", "s2 = StokesVector(1.5964694635032322,-0.2407997868272278,0.5994372371241194,-0.8510075390636394)\n", "m1 * s1 = StokesVector(1.5750484699343372,-0.46864121761363126,0.11657152710872415,0.8816680960253582)\n", "m1 @ s1 = StokesVector(1.5750484699343372,-0.46864121761363126,0.11657152710872415,0.8816680960253582)\n", "m1.dot(s1) = [ 1.57504847 -0.46864122 0.11657153 0.8816681 ]\n", "m1 @ m2 = [[ 0.88968056 0.54017069 -0.0454817 0.03081172]\n", " [-0.28616721 -0.22534394 0.15237873 0.14794223]\n", " [-0.15261965 -0.13695333 -0.13154859 0.16915209]\n", " [ 0.46629398 0.39577053 0.00381207 -0.01415688]]\n", "s1 @ m1 = StokesVector(1.062536738150666,-0.2823616179668839,0.4716926877767631,-0.5230727916610801)\n", "s1 @ s2 = 2.661695204424655\n" ] } ], "source": [ "m1 = RandomMuellerMatrix(4)\n", "m2 = RandomMuellerMatrix(4)\n", "s1 = RandomStokesVector(4)\n", "s2 = RandomStokesVector(4)\n", "print(\"m1 = \\n\",m1)\n", "print(\"m2 = \\n\",m2)\n", "print(\"s1 = \",s1)\n", "print(\"s2 = \",s2)\n", "print(\"m1 * s1 = \",m1 * s1)\n", "print(\"m1 @ s1 = \",m1 * s1)\n", "print(\"m1.dot(s1) = \",m1.dot(s1))\n", "print(\"m1 @ m2 = \",m1 * m2)\n", "print(\"s1 @ m1 = \",s1 * m1)\n", "print(\"s1 @ s2 = \",s1 * s2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mueller matrix decompositions\n", "--------------------------------\n", "\n", "So, now that we can create some physically valid Mueller matrices, we can explore some decompositions. The Cloude decomposition* decomposes a Mueller matrix into an incoherent, positive (convex) sum of Jones-Mueller matrices. \n", "\n", "*S. R. Cloude, “Group theory and polarization algebra,” Optik (Stuttgart) **75**, 26–36 (1986)\n", "and S. R. Cloude and E. Pottier, “A review of target decomposition theorems in radar polarimetry,” IEEE Trans.\n", "Geosci. and Remote Sensing 34, 498–518 (1996)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "M = \n", "\n", " [[ 0.97198146 0.30841196 -0.13263167 0.11331983]\n", " [-0.01312405 0.05122653 0.24274521 0.03582365]\n", " [-0.14679568 0.04987054 -0.08064983 -0.34434852]\n", " [-0.09716448 -0.16663341 -0.09086824 0.40486368]] \n", "\n", "Cloude decomposition =\n", "\n", "[[ 0.49815716 0.23393238 -0.05659197 0.11395561]\n", " [-0.07926307 -0.06762069 0.42302718 0.00239692]\n", " [-0.22600176 -0.30168327 -0.00377822 -0.37053662]\n", " [-0.11641932 -0.36930565 -0.03852329 0.2300661 ]] \n", "\n", "[[ 0.31877226 0.15566403 -0.09501666 0.02225629]\n", " [ 0.10148151 0.03766427 -0.2769091 0.00788727]\n", " [ 0.13742785 0.28393391 -0.01393898 -0.07703643]\n", " [ 0.06760094 0.10027555 -0.0040224 0.2497188 ]] \n", "\n", "[[ 0.15220157 -0.07995128 0.01814329 -0.02246167]\n", " [-0.03555445 0.08268129 0.09841471 0.02611275]\n", " [-0.05939849 0.06973868 -0.06428279 0.10233195]\n", " [-0.04933165 0.10311127 -0.04950611 -0.07273397]] \n", "\n", "[[ 0.00285046 -0.00123317 0.00083367 -0.0004304 ]\n", " [ 0.00021197 -0.00149834 -0.00178758 -0.00057329]\n", " [ 0.00117673 -0.00211878 0.00135016 0.00089258]\n", " [ 0.00098556 -0.00071458 0.00118356 -0.00218725]] \n", "\n", "Purity1 = 0.18455587955678526\n", "Purity2 = 0.527300463403825\n", "Purity3 = 0.9882694784032228\n" ] } ], "source": [ "M = RandomMuellerMatrix(4)\n", "MCD = M.Cloude_Decomposition()\n", "print(\"M = \\n\\n\",M,\"\\n\\nCloude decomposition =\\n\")\n", "for m in MCD:\n", " print(m,\"\\n\")\n", " \n", "print(f\"Purity1 = {(MCD[0][0,0]-MCD[1][0,0])/M[0,0]}\")\n", "print(f\"Purity2 = {(MCD[0][0,0]+MCD[1][0,0]-2*MCD[2][0,0])/M[0,0]}\")\n", "print(f\"Purity3 = {(MCD[0][0,0]+MCD[1][0,0]+MCD[2][0,0]-3*MCD[3][0,0])/M[0,0]}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that if we create the random Mueller matrix from less than four Jones-Mueller matrices, the Cloude decomposition will have only as many non-zero matrices as there were Jones-Mueller matrices. In the following, we plot the average magnitude of each element of the Cloude decomposition as a function of how many Jones-Mueller matrices were added." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-03-02T13:11:30.593692\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.5.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "jmax = 100\n", "imax = 100\n", "mCD00sums = []\n", "for i in range(1,imax):\n", " mCD00sum = np.array([0,0,0,0])\n", " for j in range(jmax):\n", " M = RandomMuellerMatrix(i)\n", " mCD = Cloude_Decomposition(M)\n", " mCD00sum = mCD00sum + np.array([m[0,0] for m in mCD])\n", " mCD00sums.append(mCD00sum/jmax)\n", "plt.figure()\n", "plt.plot(mCD00sums)\n", "plt.xlabel(\"Number of Jones-Mueller matrices averaged\")\n", "plt.ylabel(\"Cloude decomposition\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Lu-Chipman decomposition* decomposes a given Mueller matrix $\\mathbf{M}$ into the ordered product of three matrices:\n", "\n", "$$\n", "\\mathbf{M} = \\mathbf{M}_\\Delta \\mathbf{M}_R \\mathbf{M}_D\n", "$$\n", "where $\\mathbf{M}_\\Delta$ is a depolarizer, $\\mathbf{M}_R$ is a retarder, and $\\mathbf{M}_D$ is a diattenuator. Notice that if there is only one Jones-Mueller matrix used to create the random matrix, there will be no depolarization.\n", "\n", "*S.-Y. Lu and R. A. Chipman, “Interpretation of Mueller\n", "matrices based on polar decomposition,” J. Opt. Soc. Am. A \n", "13, 1106–1113 (1996)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " M =\n", "\n", "[[ 0.84442092 0.1075449 0.06352588 0.05409045]\n", " [ 0.21973306 0.44426721 -0.13529821 0.11410266]\n", " [ 0.23420748 0.03277545 -0.11947241 -0.13426277]\n", " [ 0.08729154 0.07995964 -0.03015655 -0.0153505 ]]\n", "\n", " Depolarizer =\n", "\n", "[[1. 0. 0. 0. ]\n", " [0.20185409 0.54228623 0.00714284 0.07993953]\n", " [0.30106685 0.00714284 0.24083235 0.04079128]\n", " [0.09770437 0.07993953 0.04079128 0.03708036]]\n", "\n", " Retarder =\n", "\n", "[[ 1.00000000e+00 2.10335221e-17 -1.89193279e-17 0.00000000e+00]\n", " [ 1.11022302e-16 8.54280664e-01 -3.97300271e-01 3.35197020e-01]\n", " [ 4.44089210e-16 -1.03253526e-01 -7.61690669e-01 -6.39660874e-01]\n", " [ 0.00000000e+00 5.09453881e-01 5.11839642e-01 -6.91720264e-01]]\n", "\n", " Diattenuator =\n", "\n", "[[0.84442092 0.1075449 0.06352588 0.05409045]\n", " [0.1075449 0.84027184 0.00407193 0.00346713]\n", " [0.06352588 0.00407193 0.8357836 0.00204801]\n", " [0.05409045 0.00346713 0.00204801 0.83512216]]\n", "\n", " The product Depolarizer . Retarder . Diattenuator =\n", "\n", "[[ 0.84442092 0.1075449 0.06352588 0.05409045]\n", " [ 0.21973306 0.44426721 -0.13529821 0.11410266]\n", " [ 0.23420748 0.03277545 -0.11947241 -0.13426277]\n", " [ 0.08729154 0.07995964 -0.03015655 -0.0153505 ]]\n" ] } ], "source": [ "def printMatrix(name,M):\n", " print(\"\\n\",name,\"=\\n\")\n", " print(M)\n", "\n", "M = RandomMuellerMatrix(4)\n", "depol, ret, diatten = Lu_Chipman_Decomposition(M)\n", "printMatrix(\"M\",M)\n", "printMatrix(\"Depolarizer\",depol)\n", "printMatrix(\"Retarder\",ret)\n", "printMatrix(\"Diattenuator\",diatten)\n", "printMatrix(\"The product Depolarizer . Retarder . Diattenuator\",depol @ ret @ diatten)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The reverse Lu-Chipman decomposition is similar to the Lu-Chipman decomposition, except that the order of the elements us reversed. That is,\n", "$$\n", "\\mathbf{M} = \\mathbf{M}^\\prime_D \\mathbf{M}^\\prime_R \\mathbf{M}^\\prime_\\Delta \n", "$$\n", "where $\\mathbf{M}^\\prime_D$ is a diattenuator, $\\mathbf{M}^\\prime_R$ is a retarder, and $\\mathbf{M}^\\prime_\\Delta$ is a depolarizer. Note that highly depolarizing Mueller matrices often cannot be decomposed in this fashion, an exception will be thrown if this is the case." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " M =\n", "\n", "[[ 0.79050245 0.31080475 0.13624337 -0.2012983 ]\n", " [-0.22657717 -0.64638751 -0.16880211 -0.25056315]\n", " [ 0.19461597 0.0578249 0.53578195 -0.31627667]\n", " [ 0.17512845 0.33887398 -0.24140583 -0.45722454]]\n", "\n", " Diattenuator =\n", "\n", "[[ 0.87059006 -0.230308 0.2884748 0.30995087]\n", " [-0.230308 0.75822452 -0.04163903 -0.04473893]\n", " [ 0.2884748 -0.04163903 0.77713682 0.05603824]\n", " [ 0.30995087 -0.04473893 0.05603824 0.78519151]]\n", "\n", " Retarder =\n", "\n", "[[ 1. 0. 0. 0. ]\n", " [ 0. -0.87883829 -0.22590361 -0.4202509 ]\n", " [ 0. 0.01769794 0.86476707 -0.50186123]\n", " [ 0. 0.4767914 -0.44849244 -0.75599239]]\n", "\n", " Depolarizer =\n", "\n", "[[ 1. 0. 0. 0. ]\n", " [-0.07088282 0.91237903 0.01709567 0.05714633]\n", " [-0.01846463 0.01709567 0.82188915 0.01194215]\n", " [ 0.18385422 0.05714633 0.01194215 0.79183899]]\n", "\n", " The product Diattenuator . Retarder . Depolarizer =\n", "\n", "[[ 0.79050245 0.31080475 0.13624337 -0.2012983 ]\n", " [-0.22657717 -0.64638751 -0.16880211 -0.25056315]\n", " [ 0.19461597 0.0578249 0.53578195 -0.31627667]\n", " [ 0.17512845 0.33887398 -0.24140583 -0.45722454]]\n" ] } ], "source": [ "M = RandomMuellerMatrix(2)\n", "diatten, ret, depol = Reverse_Lu_Chipman_Decomposition(M)\n", "printMatrix(\"M\",M)\n", "printMatrix(\"Diattenuator\",diatten)\n", "printMatrix(\"Retarder\",ret)\n", "printMatrix(\"Depolarizer\",depol)\n", "printMatrix(\"The product Diattenuator . Retarder . Depolarizer\", diatten @ ret @ depol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The symmetric decomposition, introduced by Ossikovski,* is a decomposition that treats the ordering of the Mueller matrix product as symmetric. That is, the Mueller matrix $\\mathbf{M}$ is decomposed into \n", "$$\n", "\\mathbf{M} = \\mathbf{M}_{D_2} \\mathbf{M}_{R_2} \\mathbf{M}_\\Delta \\mathbf{M}_{R_1} \\mathbf{M}_{D_1}\n", "$$\n", "where $\\mathbf{M}_{D_1}$ and $\\mathbf{M}_{D_2}$ are a diattenuators, $\\mathbf{M}^\\prime_{R_1}$ and $\\mathbf{M}^\\prime_{R_2}$ are retarders, and $\\mathbf{M}_\\Delta$ is a diagonal depolarizer.
\n", "\n", "*R. Ossikovski, \"Analysis of depolarizing Mueller matrices\n", "through a symmetric decomposition,\" J. Opt. Soc. Am. A, 26(5), 1109-1118 (2009)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " M =\n", "\n", "[[ 1.31346355 -0.70433204 -0.05792289 -0.42159147]\n", " [-0.59531681 0.30343934 -0.47341912 0.52095325]\n", " [ 0.31783609 -0.24582591 -0.39663839 0.37667674]\n", " [ 0.27500732 -0.26594446 0.04012768 -0.18661019]]\n", "\n", " diatten2 =\n", "\n", "[[ 1. -0.42564353 0.31261478 0.10704983]\n", " [-0.42564353 0.94073622 -0.0722223 -0.02473135]\n", " [ 0.31261478 -0.0722223 0.8954451 0.01816399]\n", " [ 0.10704983 -0.02473135 0.01816399 0.84862125]]\n", "\n", " ret2 =\n", "\n", "[[ 1. 0. 0. 0. ]\n", " [ 0. 0.70579373 -0.26909955 0.6553172 ]\n", " [ 0. 0.69389627 0.44894296 -0.56299039]\n", " [ 0. -0.14269958 0.85207726 0.5035883 ]]\n", "\n", " depol =\n", "\n", "[[1.31231491 0. 0. 0. ]\n", " [0. 1.31231491 0. 0. ]\n", " [0. 0. 0.24752227 0. ]\n", " [0. 0. 0. 0.24752227]]\n", "\n", " ret1 =\n", "\n", "[[ 1. 0. 0. 0. ]\n", " [ 0. -0.01944212 -0.73046083 0.6826778 ]\n", " [ 0. -0.9627385 0.19788894 0.18432186]\n", " [ 0. -0.26973428 -0.65365659 -0.70709015]]\n", "\n", " diatten1 =\n", "\n", "[[ 1. -0.50453294 -0.15004366 -0.31959462]\n", " [-0.50453294 0.93027935 0.04234118 0.09018718]\n", " [-0.15004366 0.04234118 0.80049588 0.02682087]\n", " [-0.31959462 0.09018718 0.02682087 0.84503274]]\n", "\n", " The product Diattenuator2 . Retarder2 . Depolarizer . Retarder1 . Diattenuator1 =\n", "\n", "[[ 1.31346355 -0.70433204 -0.05792289 -0.42159147]\n", " [-0.59531681 0.30343934 -0.47341912 0.52095325]\n", " [ 0.31783609 -0.24582591 -0.39663839 0.37667674]\n", " [ 0.27500732 -0.26594446 0.04012768 -0.18661019]]\n" ] } ], "source": [ "M = RandomMuellerMatrix(2)\n", "diatten2, ret2, depol, ret1, diatten1 = Symmetric_Decomposition(M)\n", "printMatrix(\"M\",M)\n", "printMatrix(\"diatten2\",diatten2)\n", "printMatrix(\"ret2\",ret2)\n", "printMatrix(\"depol\",depol)\n", "printMatrix(\"ret1\",ret1)\n", "printMatrix(\"diatten1\",diatten1)\n", "printMatrix(\"The product Diattenuator2 . Retarder2 . Depolarizer . Retarder1 . Diattenuator1\",\n", " diatten2 @ ret2 @ depol @ ret1 @ diatten1)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The differential decomposition is simply the matrix logarithm of the Mueller matrix." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " M =\n", "\n", "[[ 0.90254645 0.27385362 0.30832524 -0.28841948]\n", " [-0.32530387 -0.31270702 -0.36702936 -0.22067358]\n", " [-0.10609641 -0.01203455 -0.22551256 0.15574008]\n", " [-0.40125467 -0.27026781 -0.04613747 0.37130898]]\n", "\n", " L =\n", "\n", "[[-0.28858548 0.16185659 0.85000633 -0.74234137]\n", " [-1.10744528 -2.89640883 -6.81363824 0.49688726]\n", " [ 0.53841303 1.40068116 0.37418308 0.5630028 ]\n", " [-1.01346942 -0.95694393 -2.34927816 -0.91739108]]\n", "\n", " MuellerExp(L) =\n", "\n", "[[ 0.90254645 0.27385362 0.30832524 -0.28841948]\n", " [-0.32530387 -0.31270702 -0.36702936 -0.22067358]\n", " [-0.10609641 -0.01203455 -0.22551256 0.15574008]\n", " [-0.40125467 -0.27026781 -0.04613747 0.37130898]]\n" ] } ], "source": [ "M = RandomMuellerMatrix(4)\n", "L = MuellerLog(M)\n", "printMatrix(\"M\",M)\n", "printMatrix(\"L\",L)\n", "printMatrix(\"MuellerExp(L)\",MuellerExp(L))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other parameters\n", "----------------\n", "One can obtain other parameters associated with a Mueller matrix." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tmax = 1.1130487259621076\n", "Tmin = 0.3550748776407327\n", "diattenuation = 0.5162874886428319\n", "linear diattenuation = 0.5008220205047619\n", "polarization dependent loss = 4.961942308101194\n", "polarizance = 0.5633462437213892\n", "depolarization index = 0.6589112048969448\n", "extinction_ratio = 3.1346873463920426\n" ] } ], "source": [ "M = RandomMuellerMatrix(4)\n", "\n", "print(\"Tmax = \",M.Tmax())\n", "print(\"Tmin = \",M.Tmin())\n", "print(\"diattenuation = \",M.diattenuation())\n", "print(\"linear diattenuation = \",M.linear_diattenuation())\n", "print(\"polarization dependent loss = \",M.polarization_dependent_loss())\n", "print(\"polarizance = \",M.polarizance())\n", "print(\"depolarization index = \",M.depolarization_index())\n", "print(\"extinction_ratio = \",M.extinction_ratio())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is also a class `CharacterizedMueller` that performs a Lu-Chipman decomposition and extracts parameters from them. This is treated as a separate class, because initialization performs the decomposition, assigning values to the parameters. Note that the parameters differ from those shown above, due to differences in definition.\n", "\n", "N. Ghosh, M.F.G. Wood, and I.A. Vitkin, ''Mueller matrix decomposition for extraction of individual polarization parameters from complex turbid media exhibiting multiple scattering, optical activity, and linear birefringence,'' J. Biomedical Opt. 13, 044036 (2008)." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mdepol = \n", "[[ 1. 0. 0. 0. ]\n", " [-0.08476877 0.27377418 -0.03728996 0.09656993]\n", " [-0.20204637 -0.03728996 0.44373316 -0.09981802]\n", " [-0.26307903 0.09656993 -0.09981802 0.59248849]]\n", "Mret = \n", "[[ 1.00000000e+00 -9.45424294e-17 -7.28583860e-17 -5.55111512e-17]\n", " [ 0.00000000e+00 -5.62772927e-01 -5.92775720e-01 -5.76110734e-01]\n", " [-1.24900090e-16 -3.11903969e-01 7.97714968e-01 -5.16107299e-01]\n", " [ 5.55111512e-17 7.65508031e-01 -1.10759990e-01 -6.33821488e-01]]\n", "Mdiatten = \n", "[[ 0.7340618 -0.34245359 0.13371808 0.09206573]\n", " [-0.34245359 0.71472066 -0.03360347 -0.0231362 ]\n", " [ 0.13371808 -0.03360347 0.64178295 0.00903401]\n", " [ 0.09206573 -0.0231362 0.00903401 0.63488175]]\n", "DiattenuationVector = [-0.46651875 0.18216188 0.12541959]\n", "Diattenuation = 0.516287488642832\n", "CircularDiattenuation = 0.12541959458812751\n", "LinearDiattenuation = 0.5008220205047619\n", "DiattenuationAngle = 1.384664060309693 rad (79.33540669919348 deg)\n", "PolarizanceVector = [-0.08476877 -0.20204637 -0.26307903]\n", "Polarizance = 0.342372684079807\n", "DepolarizationCoefficient = 0.4366652735461514\n", "LinearRetardance = 2.2572802432682404 rad (129.33263111753394 deg)\n", "OpticalRotation = 0.4371035715152117 rad (25.04418985791638 deg)\n", "Retardance = 2.345409580619083 rad (134.38207019902185deg)\n", "RetardanceAngle = 0.7135526418560333 rad (40.88355483876068 deg)\n" ] } ], "source": [ "parametersM = CharacterizedMueller(M)\n", "print(parametersM)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.13" } }, "nbformat": 4, "nbformat_minor": 2 }