Testing, and some streamlining, of LUT calibration (& CLI)

hrl/calibration/measurement.py

@@ -1,4 +1,5 @@
 from functools import partial
+from pathlib import Path
 
 import numpy as np
 
@@ -55,12 +56,14 @@ def draw_uniform_square(ihrl, intensity, patch_size=0.5):
     )
     patch = ihrl.graphics.newTexture(np.array([[intensity]]))
     patch.draw(patch_position, (patch_width, patch_height))
+    ihrl.graphics.flip()
 
 
 def measure_lut(
     ihrl,
     intensities=setup_intensities(0.0, 1.0, 2**16),
     stim_draw_func=partial(draw_uniform_square, patch_size=0.5),
+    out_file=None,
     n_samples=5,
     sleep_time=200,
 ):
@@ -75,35 +78,47 @@ def measure_lut(
     stim_draw_func : callable, optional
         function with signature `(ihrl, intensity)` that draws the stimulus
         for each measurement; defaults to `draw_uniform_square(patch_size=0.5)`
+    out_file : str or Path, optional
+        path to output file for measurements, by default None (no file output)
     n_samples : int
         number of photometer readings per intensity level, by default 5
     sleep_time : float
         time (ms) to wait between photometer readings, by default 200ms
     """
+    measurements = np.ndarray((len(intensities), n_samples + 1), dtype=float)
 
     for idx_int, intensity in enumerate(intensities):
         print(
             f"Current Intensity: {intensity:.2f} "
             f"[{idx_int:d} of {len(intensities)} "
             f"({idx_int / len(intensities) * 100:.1f}%)]"
         )
-        ihrl.results["Intensity"] = intensity
 
         # Draw (update) stimulus
         stim_draw_func(ihrl, intensity)
-        ihrl.graphics.flip()
 
         # Multiple samples for each intensity value
         for idx_sample in range(n_samples):
             sample = ihrl.photometer.readLuminance(5, int(sleep_time))
-            ihrl.results[f"Luminance{idx_sample}"] = sample
+            measurements[idx_int, idx_sample + 1] = sample
 
         # Write measured samples to file
-        ihrl.writeResultLine()
-
-        if ihrl.inputs.checkEscape():
+        measurements[idx_int, 0] = intensity
+        if out_file is not None:
+            out_file = Path(out_file).expanduser().resolve()
+            np.savetxt(
+                out_file,
+                measurements,
+                delimiter=",",
+                header="intensity," + ",".join([f"luminance{i}" for i in range(n_samples)]),
+                comments="",
+            )
+
+        if ihrl.inputs is not None and ihrl.inputs.checkEscape():
             break
 
+    return measurements
+
 
 def combine(measurements):
     """Construct an intensity-to-luminance map from (sets of) measurements

hrl/luts.py

@@ -125,12 +125,16 @@ def create_lut(
             intensity_out = intensity_in^(1/gamma)
             luminance = k * intensity_in^gamma + dark
     """
+    # Linearly spaced input intensities from 0 to 1
     x = np.linspace(0.0, 1.0, n)
 
-    out = x ** (1 / gamma)
+    # Linearly spaced luminance value.
+    # Use linspace, because that's also what we use in the calibration (linearize())
+    # This way, luminance values are bit-identical, important for testing and reproducibility.
+    lum = np.linspace(dark, k + dark, n)
 
-    lum = k * x**gamma + dark
-    lum[0] = dark  # enforce zero row
+    # Output intensities: apply inverse gamma correction
+    out = x ** (1 / gamma)
 
     return np.column_stack([x, out, lum])
 

hrl/photometer/mock.py

@@ -0,0 +1,97 @@
+"""Mock photometer (for testing) using a lookup table to simulate luminance readings.
+
+Simulated photometer object that adheres to the same interface as the real photometer,
+but returns simulated luminance values (based on a LUT) instead of reading from a physical device.
+
+The challenge with mocking a photometer is that it has no knowledge of what
+intensity is currently being displayed -- that information travels through the
+physical path (screen → photons → sensor).
+
+MockPhotometer solves this by holding ``current_intensity`` as a direct attribute,
+which a stimulus draw function should update before ``readLuminance`` is called.
+
+Typical usage::
+
+    ihrl.photometer = MockPhotometer(lut)
+
+    def mock_draw(ihrl, intensity):
+        ihrl.photometer.current_intensity = intensity
+        # draw_uniform_square(ihrl, intensity)
+
+    measure_lut(ihrl, stim_draw_func=mock_draw, ...)
+"""
+
+import numpy as np
+
+from .photometer import Photometer
+
+
+class MockPhotometer(Photometer):
+    """Photometer that returns simulated luminance values from a LUT.
+
+    The currently displayed intensity is tracked via the ``current_intensity``
+    attribute, which must be updated by the stimulus drawing code before each
+    call to ``readLuminance``.
+
+    Parameters
+    ----------
+    lut : array-like or callable
+        Defines the intensity to luminance mapping.
+
+        - **array**: shape ``(N, 2)`` or ``(N, 3)``.  Column 0 is intensity
+          (input), the *last* column is luminance (cd/m²).  Values are
+          linearly interpolated.
+        - **callable**: ``lut(intensity: float) -> float`` returning cd/m².
+    noise : float, optional
+        Standard deviation of zero-mean Gaussian noise added to each reading,
+        in cd/m², by default 0.0 (noiseless).
+    rng : numpy.random.Generator or int or None, optional
+        Random number generator (or seed) used for noise.  Pass an integer
+        for reproducible results. By default None (unseeded).
+
+    Attributes
+    ----------
+    current_intensity : float
+        The intensity value (in [0, 1]) most recently drawn to the screen.
+        Update this before calling ``readLuminance``.
+    """
+
+    def __init__(self, lut, noise=0.0, rng=None):
+        super().__init__()
+        self.current_intensity = 0.0
+        self.noise = noise
+        self.rng = np.random.default_rng(rng)
+
+        if callable(lut):
+            self._lut_func = lut
+        else:
+            lut = np.asarray(lut)
+            if lut.shape[1] == 3:
+                # Full 3-column LUT (intensity_in, intensity_out, luminance):
+                # use intensity_out (col 1) as the physical intensity axis.
+                intensities = lut[:, 1]
+            else:
+                intensities = lut[:, 0]
+            luminances = lut[:, -1]
+            self._lut_func = lambda x: float(np.interp(x, intensities, luminances))
+
+    def readLuminance(self, n=3, slp=None):
+        """Return simulated luminance for the currently displayed intensity.
+
+        Parameters
+        ----------
+        n : int
+            Number of samples to average (mirrors the real photometer API).
+        slp : int
+            Sleep time between samples in ms (ignored in mock).
+
+        Returns
+        -------
+        float
+            Simulated luminance in cd/m², averaged over ``n`` samples.
+        """
+        base_lum = self._lut_func(self.current_intensity)
+        if self.noise > 0.0:
+            samples = base_lum + self.rng.normal(0.0, self.noise, size=n)
+            return float(np.mean(samples))
+        return base_lum

hrl/util/lut/linearize.py

@@ -1,4 +1,5 @@
 import argparse
+from pathlib import Path
 
 import numpy as np
 
@@ -18,25 +19,40 @@
     add_help=False,
     parents=[intensities_argparser],
 )
+parser.add_argument(
+    "-i",
+    "--in_file",
+    default="smooth.csv",
+    type=Path,
+    help="path to input measurements csv, by default 'smooth.csv'",
+)
+parser.add_argument(
+    "-o",
+    "--out_file",
+    default="lut.csv",
+    type=Path,
+    help="path for output csv, by default 'lut.csv'",
+)
 
 
 def command(parsed_args):
     """Sample a linear subset of the gamma table"""
-    n_steps = 2**parsed_args.bit_depth
 
     # Load (smoothed) LUT
-    lut = np.genfromtxt("smooth.csv", skip_header=1, delimiter=",")
+    in_file = parsed_args.in_file.expanduser().resolve()
+    print(f"Loading measurements from {in_file} ...")
+    measurements = np.genfromtxt(in_file, delimiter=",", skip_header=1)
 
     # Linearize LUT
-    linearized_lut = hrl.calibration.measurement.linearize(lut, bit_depth=parsed_args.bit_depth)
-
-    # Save linearized LUT to file
-    print("Saving to File...")
-    out_file = open("lut.csv", "w")
-    headers = "intensity_in,intensity_out,luminance\n"
-    out_file.write(headers)
-    np.savetxt(out_file, linearized_lut, delimiter=",")
-    out_file.close()
+    linearized_lut = hrl.calibration.measurement.linearize(
+        measurements, bit_depth=parsed_args.bit_depth
+    )
+
+    # Write to file
+    out_file = parsed_args.out_file.expanduser().resolve()
+    print(f"Saving to {out_file} ...")
+    headers = "intensity_in,intensity_out,luminance"
+    np.savetxt(out_file, linearized_lut, delimiter=",", header=headers, comments="")
 
 
 if __name__ == "__main__":

hrl/util/lut/measure.py

@@ -1,6 +1,7 @@
 import argparse
 from datetime import timedelta
 from functools import partial
+from pathlib import Path
 from timeit import default_timer as timer
 
 from hrl import HRL
@@ -58,9 +59,9 @@
 parser.add_argument(
     "-o",
     "--out_file",
-    type=str,
+    type=Path,
     default="measure.csv",
-    help="Output filename, by default 'measure.csv'",
+    help="path to output measurements csv, by default 'measure.csv'",
 )
 
 
@@ -71,8 +72,6 @@ def command(parsed_args):
     start = timer()
 
     # Initializing HRL
-    headers = ["Intensity"] + ["Luminance" + str(i) for i in range(parsed_args.n_samples)]
-
     ihrl = HRL(
         graphics=parsed_args.graphics,
         inputs="keyboard",
@@ -84,8 +83,6 @@ def command(parsed_args):
         wdth_offset=parsed_args.width_offset,
         db=True,
         scrn=parsed_args.screen,
-        rfl=parsed_args.out_file,
-        rhds=headers,
     )
 
     # Set up intensity values to be measured
@@ -105,6 +102,7 @@ def command(parsed_args):
         ihrl,
         intensities=intensities,
         stim_draw_func=partial(draw_uniform_square, patch_size=parsed_args.patch_size),
+        out_file=parsed_args.out_file,
         n_samples=parsed_args.n_samples,
         sleep_time=parsed_args.sleep_time,
     )

hrl/util/lut/plot.py

@@ -1,4 +1,5 @@
 import argparse
+from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -11,10 +12,18 @@
     """,
     add_help=False,
 )
+parser.add_argument(
+    "-i",
+    "--in_file",
+    default="lut.csv",
+    type=Path,
+    help="path to LUT csv plot, by default 'lut.csv'",
+)
 
 
 def command(parsed_args):
-    lut = np.genfromtxt("lut.csv", skip_header=1, delimiter=",")
+    lut_file = parsed_args.in_file.expanduser().resolve()
+    lut = np.genfromtxt(lut_file, delimiter=",", skip_header=1)
 
     plt.figure()
 

hrl/util/lut/smooth.py

@@ -1,4 +1,5 @@
 import argparse
+from pathlib import Path
 
 import numpy as np
 
@@ -36,12 +37,33 @@
     nargs="+",
     help="kernel for smoothing, by default [0.2, 0.2, 0.2, 0.2, 0.2]",
 )
+parser.add_argument(
+    "-i",
+    "--in_file",
+    default="measure.csv",
+    type=Path,
+    nargs="+",
+    help="path(s) for input measurement csv(s), by default 'measure.csv'",
+)
+parser.add_argument(
+    "-o",
+    "--out_file",
+    default="smooth.csv",
+    type=Path,
+    help="path for output smoothed measurements csv, by default 'smooth.csv'",
+)
 
 
 def command(parsed_args):
     # Load measurement data
-    files = ["measure.csv"]
-    measurements = [np.genfromtxt(fl, skip_header=1, delimiter=",") for fl in files]
+    measurements = []
+    in_files = parsed_args.in_file
+    if isinstance(in_files, (str, Path)):
+        in_files = [in_files]
+    for file in in_files:
+        filename = file.expanduser().resolve()
+        print(f"Loading from {filename} ...")
+        measurements.append(np.genfromtxt(filename, delimiter=",", skip_header=1))
 
     # Combine
     luminance_map = hrl.calibration.measurement.combine(measurements)
@@ -58,12 +80,10 @@ def command(parsed_args):
     )
 
     # Save smoothed LUT to file
-    print("Saving to File...")
-    out_filename = "smooth.csv"
-    out_file = open(out_filename, "w")
-    out_file.write("intensity_in,luminance\n")
-    np.savetxt(out_file, table, delimiter=",")
-    out_file.close()
+    out_file = parsed_args.out_file.expanduser().resolve()
+    print(f"Saving to {out_file}...")
+    header = "intensity_in,luminance"
+    np.savetxt(out_file, table, delimiter=",", header=header, comments="")
 
 
 if __name__ == "__main__":