Add intensity gradient compass.

Author: evgkanias

templates/plot_all.py

@@ -11,8 +11,8 @@
 import matplotlib.pyplot as plt
 
 if __name__ == "__main__":
-    data_dir = os.path.join(os.getcwd(),'../sardinia_data')
-    out_base_dir = os.path.join(os.getcwd(),'../plots')
+    data_dir = os.path.join(os.getcwd(), '../sardinia_data')
+    out_base_dir = os.path.join(os.getcwd(), '../plots')
 
     # Make plot directory if it doesn't exist. Existing files will
     # be overwritten.

templates/plot_from_dict.py

@@ -16,6 +16,7 @@
 import os
 from scipy.stats import circmean
 
+
 def act(s, log=True):
     """
     Photoreceptor activation function. Gkanias et al. (2019) uses square
@@ -31,6 +32,7 @@ def act(s, log=True):
 
     return np.log(s)
 
+
 def synchronise(d, zero=True, unwrap=True):
     """
     Zero each angle series and adjust values so wraps plot nicely.
@@ -53,10 +55,9 @@ def synchronise(d, zero=True, unwrap=True):
 
     return d
 
-def decode_sky_compass(po,
-                       n_sol=8,
-                       pol_prefs=np.radians([0,90,180,270,45,135,225,315])
-):
+
+def decode_sky_compass(po, n_sol=8, pol_prefs=np.radians([0, 90, 180, 270, 45, 135, 225, 315]),
+                       polarisation=True, intensity=False):
     """
     Sky compass decoding routine from Gkanias et al. (2019). This has been
     implemeneted from scratch using the equations presented as a demonstration
@@ -77,65 +78,93 @@ def decode_sky_compass(po,
     # po[p][t][d] reading for diode d of unit p at time t
     units, duration, diodes = po.shape
 
-    n_pol = len(pol_prefs)
-    sol_prefs = np.linspace(0,2*np.pi - (2*np.pi / n_sol),n_sol)
+    sol_prefs = np.linspace(0, 2 * np.pi - (2 * np.pi / n_sol), n_sol)
 
     angular_outputs = []
     confidence_outputs = []
 
-    for t in range(duration): # For each timestep
-        # Get photodiode responses for each unit
-        unit_responses = [ po[x][t] for x in range(units) ]
-        ss = [ (s_h, s_v) for [_, _, s_v, s_h] in unit_responses ]
-        rates = [ (act(h), act(v)) for (h, v) in ss ]
-        r_pol = [ (r_h - r_v) / (r_h + r_v) for (r_h, r_v) in rates ]
-
-        # Init for sum
-        r_sol = np.zeros(n_sol)
-        R = 0
-        for z in range(n_sol):
-            for j in range(n_pol): # Compute SOL responses (Eq. (2))
-                aj = pol_prefs[j] - np.pi/2
-                r_sol[z] += (n_sol/n_pol) * np.sin(aj - sol_prefs[z]) * r_pol[j]
-
-            # Accumulate R (Eq. (3))
-            R += r_sol[z]*np.exp(-1j*2*np.pi*(z - 1) / n_sol)
-
-        a = np.real(R)
-        b = np.imag(R)
-
-        # Compute argument and magnitude of complex conjugate of R (Eq. (4))
-        angle = np.arctan2(-b,a)
-        confidence = np.sqrt(a**2 + b**2)
+    for t in range(duration):  # For each timestep
+
+        response = np.clip(po[:, t] / 12000., 0., 1.)
+        if polarisation and not intensity:
+            angle, sigma = pol2sol(response, sol_prefs, pol_prefs, unit_tranform=unit2pol)
+        elif not polarisation and intensity:
+            angle, sigma = pol2sol(response, sol_prefs + np.pi, pol_prefs, unit_tranform=unit2int)
+        else:
+            a_pol, c_pol = pol2sol(response, sol_prefs, pol_prefs, unit_tranform=unit2pol)
+            a_int, c_int = pol2sol(response, sol_prefs + np.pi, pol_prefs, unit_tranform=unit2int)
+            s_pol = 1. / c_pol ** 2
+            s_int = 1. / c_int ** 2
+
+            angle = (s_pol * a_int + s_int * a_pol) / (s_pol + s_int)
+            sigma = (s_pol * s_int) / (s_pol + s_int)
         angular_outputs.append(angle)
-        confidence_outputs.append(confidence)
+        confidence_outputs.append(sigma)
 
     return angular_outputs, confidence_outputs
 
 
+def pol2sol(po, sol_prefs, pol_prefs, unit_tranform):
+    units, diodes = po.shape
+    n_sol = len(sol_prefs)
+    n_pol = len(pol_prefs)
+
+    # Get photodiode responses for each unit
+    r_pol = np.array([unit_tranform(po[x]) for x in range(units)])
+
+    # Init for sum
+    r_sol = np.zeros(n_sol)
+    R = 0
+    for z in range(n_sol):
+        for j in range(n_pol):  # Compute SOL responses (Eq. (2))
+            aj = pol_prefs[j] - np.pi / 2
+            r_sol[z] += (n_sol / n_pol) * np.sin(aj - sol_prefs[z]) * r_pol[j]
+
+        # Accumulate R (Eq. (3))
+        R += r_sol[z] * np.exp(-1j * 2 * np.pi * (z - 1) / n_sol)
+
+    a = np.real(R)
+    b = np.imag(R)
+
+    # Compute argument and magnitude of complex conjugate of R (Eq. (4))
+    angle = np.arctan2(-b, a)
+    confidence = np.sqrt(a ** 2 + b ** 2)
+
+    return angle, confidence
+
+
+def unit2pol(unit_response, log=True):
+    r_h, r_v = unit2vh(unit_response, log)
+    return (r_h - r_v) / (r_h + r_v)
+
+
+def unit2int(unit_response, log=True):
+    r_h, r_v = unit2vh(unit_response, log)
+    r = (r_h + r_v) / 2
+    return r
+
+
+def unit2vh(unit_response, log=True):
+    _, _, s_v, s_h = unit_response
+    return act(s_h, log), act(s_v, log)
+
+
 def produce_plot(data_dictionary):
 
-    mosaic = [["image","skycompass"],
-              ["image","imu_yaw"],
-              ["image","averages"]]
+    mosaic = [["image", "skycompass"],
+              ["image", "imu_yaw"],
+              ["image", "averages"]]
 
-    fig, ax = plt.subplot_mosaic(mosaic, figsize=(10,4))
+    fig, ax = plt.subplot_mosaic(mosaic, figsize=(10, 4))
 
     #
     # Plot image if available
     #
-    imagefile = data_dictionary["image_filename"]
-    del data_dictionary["image_filename"]
-    if imagefile != None:
-        img = mpimg.imread(imagefile)
-        ax["image"].imshow(img)
-        ax["image"].set_title(imagefile.split(".")[0].replace("_", " "))
-    else:
-        ax["image"].set_title("No image found")
+    plot_image(data_dictionary["image_filename"], ax=ax["image"])
 
-    max_length = 0 # Maximum recording length in set
-    min_length = 0 # Minimum recording length in set
-    first = True # First iteration flag
+    max_length = 0  # Maximum recording length in set
+    min_length = 0  # Minimum recording length in set
+    first = True  # First iteration flag
 
     # Corrected datastructures for computing averages
     corrected_skycompass_table = []
@@ -152,7 +181,7 @@ def produce_plot(data_dictionary):
             po.append(data[k])
 
         po = np.array(po)
-        pol_sensor_angle = synchronise(np.array(decode_sky_compass(po)[0]))
+        pol_sensor_angle = synchronise(np.array(decode_sky_compass(po, polarisation=True, intensity=True)[0]))
 
         if max_length < len(imu):
             max_length = len(imu)
@@ -212,3 +241,16 @@ def produce_plot(data_dictionary):
 
     return fig
 
+
+def plot_image(imagefile, ax=None):
+
+    if ax is None:
+        ax = plt.subplot(111)
+
+    if imagefile is not None:
+        img = mpimg.imread(imagefile)
+        ax.imshow(img)
+        ax.set_title(imagefile.split(".")[0].replace("_", " "))
+    else:
+        ax.set_title("No image found")
+

templates/plot_json.py

@@ -15,8 +15,9 @@
 
 import plot_from_dict as pfd
 
+
 def plot_json(session, outfile=None):
-        """
+    """
     Produce a figure along with a local output filename. For the output filename
     the caller must specify any preceeding directory structure, otherwise the
     output file will be relative to the calling directory.
@@ -33,6 +34,7 @@ def plot_json(session, outfile=None):
     #
     # Pull in data
     #
+
     if not os.path.isdir(session):
         print("-s must specify a directory")
         sys.exit()
@@ -46,7 +48,6 @@ def plot_json(session, outfile=None):
         print("There are no json files in the specified directory.")
         sys.exit()
 
-
     if len(imagefile) == 0:
         print("Warning: missing sky image file from the session directory.")
         imagefile = None
@@ -56,14 +57,14 @@ def plot_json(session, outfile=None):
         imagefile = imagefile[0]
         # If outfile not specified, set to session spec
         if outfile == None:
-            outfile = imagefile.split(".")[0] + ".pdf"
+            outfile = imagefile.split(".")[0] + ".png"
 
     # No output file specified and no/multiple image files
     # Use session directory name as pdf name
     if outfile == None:
         session_path = session.split("/")
         print(session_path)
-        outfile = session_path[len(session_path) - 2] + ".pdf"
+        outfile = session_path[len(session_path) - 2] + ".png"
 
     # Pull all recordings into dictionary
     full_data = dict()
@@ -78,7 +79,8 @@ def plot_json(session, outfile=None):
     fig = pfd.produce_plot(full_data)
     return fig, outfile
 
-if __name__=="__main__":
+
+if __name__ == "__main__":
     calling_directory = os.getcwd()
     parser = argparse.ArgumentParser(
         description="Produce a recording-session plot from json."