filmsim/filmcolor

#!/usr/bin/env -S uv run --script
# /// script
# dependencies = [
#   "numpy",
#   "scipy",
#   "Pillow",
#   "imageio",
#   "rawpy",
#   "colour-science",
# ]
# ///

# -*- coding: utf-8 -*-
"""
Single-file Python script for Film Stock Color Emulation based on Datasheet CSV.

Focuses on color transformation, applying effects derived from datasheet parameters.
Assumes input image is in linear RGB format. Excludes film grain simulation.

Dependencies: numpy, imageio, scipy
Installation: pip install numpy imageio scipy Pillow
"""

import argparse
import csv
import numpy as np
import imageio.v3 as iio
from scipy.interpolate import interp1d
from scipy.ndimage import gaussian_filter, gaussian_filter1d
import rawpy
import os
import sys
import colour
from colour.colorimetry import SDS_ILLUMINANTS
import json

# --- Configuration ---
# Small epsilon to prevent log(0) or division by zero errors
EPSILON = 1e-10

from typing import Optional, List


class Info:
    name: str
    description: str
    format_mm: int
    version: str

    def __init__(
        self, name: str, description: str, format_mm: int, version: str
    ) -> None:
        self.name = name
        self.description = description
        self.format_mm = format_mm
        self.version = version

    def __repr__(self) -> str:
        return (
            f"Info(name={self.name}, description={self.description}, "
            f"format_mm={self.format_mm}, version={self.version})"
        )


class Balance:
    r_shift: float
    g_shift: float
    b_shift: float

    def __init__(self, r_shift: float, g_shift: float, b_shift: float) -> None:
        self.r_shift = r_shift
        self.g_shift = g_shift
        self.b_shift = b_shift
    
    def __repr__(self) -> str:
        return (
            f"Balance(r_shift={self.r_shift:.3f}, g_shift={self.g_shift:.3f}, b_shift={self.b_shift:.3f})"
        )


class Gamma:
    r_factor: float
    g_factor: float
    b_factor: float

    def __init__(self, r_factor: float, g_factor: float, b_factor: float) -> None:
        self.r_factor = r_factor
        self.g_factor = g_factor
        self.b_factor = b_factor

    def __repr__(self) -> str:
        return (
            f"Gamma(r_factor={self.r_factor:.3f}, g_factor={self.g_factor:.3f}, b_factor={self.b_factor:.3f})"
        )


class Processing:
    gamma: Gamma
    balance: Balance

    def __init__(self, gamma: Gamma, balance: Balance) -> None:
        self.gamma = gamma
        self.balance = balance

    def __repr__(self) -> str:
        return (
            f"Processing(gamma=({self.gamma.r_factor:.3f}, {self.gamma.g_factor:.3f}, {self.gamma.b_factor:.3f}), "
            f"balance=({self.balance.r_shift:.3f}, {self.balance.g_shift:.3f}, {self.balance.b_shift:.3f}))"
        )


class Couplers:
    saturation_amount: float
    dir_amount_rgb: List[float]
    dir_diffusion_um: float
    dir_diffusion_interlayer: float

    def __init__(self, saturation_amount: float, dir_amount_rgb: List[float], dir_diffusion_um: float, dir_diffusion_interlayer: float) -> None:
        self.saturation_amount = saturation_amount
        self.dir_amount_rgb = dir_amount_rgb
        self.dir_diffusion_um = dir_diffusion_um
        self.dir_diffusion_interlayer = dir_diffusion_interlayer


class HDCurvePoint:
    d: Optional[float]
    r: float
    g: float
    b: float

    def __init__(self, d: Optional[float], r: float, g: float, b: float) -> None:
        self.d = d
        self.r = r
        self.g = g
        self.b = b

    def __repr__(self) -> str:
        return f"HDCurvePoint(d={self.d}, r={self.r:.3f}, g={self.g:.3f}, b={self.b:.3f})"
    

class SpectralSensitivityCurvePoint:
    wavelength: float
    y: float
    m: float
    c: float

    def __init__(self, wavelength: float, y: float, m: float, c: float) -> None:
        self.wavelength = wavelength
        self.y = y
        self.m = m
        self.c = c

    def __repr__(self) -> str:
        return f"SpectralSensitivityCurvePoint(wavelength={self.wavelength:.1f}, y={self.y:.3f}, m={self.m:.3f}, c={self.c:.3f})"


class RGBValue:
    r: float
    g: float
    b: float

    def __init__(self, r: float, g: float, b: float) -> None:
        self.r = r
        self.g = g
        self.b = b

    def __repr__(self) -> str:
        return f"RGBValue(r={self.r:.3f}, g={self.g:.3f}, b={self.b:.3f})"


class Curves:
    hd: List[HDCurvePoint]
    spectral_sensitivity: List[SpectralSensitivityCurvePoint]

    def __init__(self, hd: List[HDCurvePoint], spectral_sensitivity: List[SpectralSensitivityCurvePoint]) -> None:
        self.hd = hd
        self.spectral_sensitivity = spectral_sensitivity

    def __repr__(self) -> str:
        return f"Curves(hd={',\n'.join(repr(point) for point in self.hd)}, spectral_sensitivity={',\n'.join(repr(point) for point in self.spectral_sensitivity)})"


class Halation:
    strength: RGBValue
    size_um: RGBValue

    def __init__(self, strength: RGBValue, size_um: RGBValue) -> None:
        self.strength = strength
        self.size_um = size_um

    def __repr__(self) -> str:
        return f"Halation(strength={self.strength}, size_um={self.size_um})"


class Interlayer:
    diffusion_um: float

    def __init__(self, diffusion_um: float) -> None:
        self.diffusion_um = diffusion_um

    def __repr__(self) -> str:
        return f"Interlayer(diffusion_um={self.diffusion_um:.1f})"


class Calibration:
    iso: int
    middle_gray_logE: float

    def __init__(self, iso: int, middle_gray_logE: float) -> None:
        self.iso = iso
        self.middle_gray_logE = middle_gray_logE

    def __repr__(self) -> str:
        return (
            f"Calibration(iso={self.iso}\nmiddle_gray_logE={self.middle_gray_logE:.3f})"
        )


class Properties:
    halation: Halation
    couplers: Couplers
    interlayer: Interlayer
    curves: Curves
    calibration: Calibration

    def __init__(
        self,
        halation: Halation,
        couplers: Couplers,
        interlayer: Interlayer,
        curves: Curves,
        calibration: Calibration,
    ) -> None:
        self.halation = halation
        self.couplers = couplers
        self.interlayer = interlayer
        self.curves = curves
        self.calibration = calibration

    def __repr__(self) -> str:
        return (
            f"Properties(halation={self.halation}\ncouplers={self.couplers}\n"
            f"interlayer={self.interlayer}\ncurves={self.curves}\n"
            f"calibration={self.calibration})"
        )


class FilmDatasheet:
    info: Info
    processing: Processing
    properties: Properties

    def __init__(
        self, info: Info, processing: Processing, properties: Properties
    ) -> None:
        self.info = info
        self.processing = processing
        self.properties = properties

    def __repr__(self) -> str:
        return (
            f"FilmDatasheet(info={self.info}\nprocessing={self.processing}\n"
            f"properties={self.properties})"
        )

import pprint

def parse_datasheet_json(json_filepath) -> FilmDatasheet | None:
    # Parse JSON into FilmDatasheet object
    """Parses the film datasheet JSON file.
    Args:
        json_filepath (str): Path to the datasheet JSON file.
    Returns:
        FilmDatasheet: Parsed datasheet object.
    """
    if not os.path.exists(json_filepath):
        print(f"Error: Datasheet file not found at {json_filepath}", file=sys.stderr)
        return None
    try:
        with open(json_filepath, "r") as jsonfile:
            data = json.load(jsonfile)
            # Parse the JSON data into the FilmDatasheet structure
            info = Info(
                name=data["info"]["name"],
                description=data["info"]["description"],
                format_mm=data["info"]["format_mm"],
                version=data["info"]["version"],
            )
            gamma = Gamma(
                r_factor=data["processing"]["gamma"]["r_factor"],
                g_factor=data["processing"]["gamma"]["g_factor"],
                b_factor=data["processing"]["gamma"]["b_factor"],
            )
            balance = Balance(
                r_shift=data["processing"]["balance"]["r_shift"],
                g_shift=data["processing"]["balance"]["g_shift"],
                b_shift=data["processing"]["balance"]["b_shift"],
            )
            processing = Processing(gamma=gamma, balance=balance)
            halation = Halation(
                strength=RGBValue(
                    r=data["properties"]["halation"]["strength"]["r"],
                    g=data["properties"]["halation"]["strength"]["g"],
                    b=data["properties"]["halation"]["strength"]["b"],
                ),
                size_um=RGBValue(
                    r=data["properties"]["halation"]["size_um"]["r"],
                    g=data["properties"]["halation"]["size_um"]["g"],
                    b=data["properties"]["halation"]["size_um"]["b"],
                ),
            )
            couplers = Couplers(
                saturation_amount=data["properties"]["couplers"]["saturation_amount"],
                dir_amount_rgb=data["properties"]["couplers"]["dir_amount_rgb"],
                dir_diffusion_um=data["properties"]["couplers"]["dir_diffusion_um"],
                dir_diffusion_interlayer=data["properties"]["couplers"]["dir_diffusion_interlayer"],
            )
            interlayer = Interlayer(
                diffusion_um=data["properties"]["interlayer"]["diffusion_um"]
            )
            calibration = Calibration(
                iso=data["properties"]["calibration"]["iso"],
                middle_gray_logE=data["properties"]["calibration"]["middle_gray_logh"],
            )
            curves = Curves(
                hd=[
                    HDCurvePoint(d=point["d"], r=point["r"], g=point["g"], b=point["b"])
                    for point in data["properties"]["curves"]["hd"]
                ],
                spectral_sensitivity=[
                    SpectralSensitivityCurvePoint(
                        wavelength=point["wavelength"],
                        y=point["y"],
                        m=point["m"],
                        c=point["c"],
                    )
                    for point in data["properties"]["curves"]["spectral_sensitivity"]
                ],
            )
            print(f"Parsed {len(curves.hd)} H&D curve points.")
            properties = Properties(
                calibration=calibration,
                halation=halation,
                couplers=couplers,
                interlayer=interlayer,
                curves=curves,
            )
            return FilmDatasheet(
                info=info, processing=processing, properties=properties
            )
    except Exception as e:
        print(f"Error parsing datasheet JSON '{json_filepath}': {e}", file=sys.stderr)
        return None


def um_to_pixels(sigma_um, image_width_px, film_format_mm):
    """Converts sigma from micrometers to pixels."""
    if film_format_mm <= 0 or image_width_px <= 0:
        return 0
    microns_per_pixel = (film_format_mm * 1000.0) / image_width_px
    sigma_pixels = sigma_um / microns_per_pixel
    return sigma_pixels



def compute_inhibitor_matrix(amount_rgb: List[float], diffusion_interlayer: float) -> np.ndarray:
    """
    Computes the 3x3 DIR coupler inhibitor matrix.
    Diagonal elements represent self-inhibition (intra-layer).
    Off-diagonal elements represent cross-inhibition (inter-layer).
    """
    # Start with an identity matrix representing the source of inhibitors
    matrix = np.eye(3)
    # Apply 1D blur across the layer axis to simulate diffusion between layers
    if diffusion_interlayer > 0:
        matrix = gaussian_filter1d(matrix, sigma=diffusion_interlayer, axis=0, mode='constant', cval=0)
        matrix = gaussian_filter1d(matrix, sigma=diffusion_interlayer, axis=1, mode='constant', cval=0)
    # Normalize rows to ensure diffusion doesn't create/destroy inhibitor
    row_sums = matrix.sum(axis=1)
    matrix = matrix / row_sums[:, np.newaxis]
    # Scale by the amount of inhibitor released by each source layer
    matrix = matrix * np.array(amount_rgb)[:, np.newaxis]
    return matrix


def compute_uncoupled_hd_curves(hd_curve_data: List[HDCurvePoint], inhibitor_matrix: np.ndarray) -> List[HDCurvePoint]:
    """
    Pre-calculates a new set of H&D curves that represent the film's
    response *without* DIR couplers.
    """
    log_E_values = np.array([p.d for p in hd_curve_data if p.d is not None])
    density_r = np.array([p.r for p in hd_curve_data])
    density_g = np.array([p.g for p in hd_curve_data])
    density_b = np.array([p.b for p in hd_curve_data])

    # For a neutral gray ramp, we assume the density forming in all three layers can be
    # approximated by the green channel's response. This is our source of inhibitors.
    neutral_density_curve = density_g

    # Create a (num_points, 3) source density matrix for the neutral ramp.
    # We tile the single neutral curve across all three source channels.
    source_density_matrix = np.tile(neutral_density_curve[:, np.newaxis], (1, 3))

    # Calculate inhibitor signal received by each destination layer.
    # This is the matrix product of the source densities and the inhibitor matrix.
    # The (i, j)-th element of inhibitor_matrix is the effect FROM source i ON destination j.
    # Shape: (40, 3) @ (3, 3) -> (40, 3)
    inhibitor_effect = source_density_matrix @ inhibitor_matrix

    # For each channel, find the new "uncoupled" density curve by shifting the
    # exposure axis by the calculated inhibitor effect for that channel.
    uncoupled_r = np.interp(log_E_values, log_E_values - inhibitor_effect[:, 0], density_r)
    uncoupled_g = np.interp(log_E_values, log_E_values - inhibitor_effect[:, 1], density_g)
    uncoupled_b = np.interp(log_E_values, log_E_values - inhibitor_effect[:, 2], density_b)

    # Reassemble into the HDCurvePoint list structure
    uncoupled_curve = []
    for i, log_e in enumerate(log_E_values):
        uncoupled_curve.append(HDCurvePoint(d=log_e, r=uncoupled_r[i], g=uncoupled_g[i], b=uncoupled_b[i]))
        
    return uncoupled_curve


def apply_dir_coupler_simulation(log_exposure_rgb, naive_density_rgb, inhibitor_matrix, diffusion_um, film_format_mm, image_width_px):
    """
    Applies the DIR coupler effect to the log exposure image.
    """
    # 1. Spatially diffuse the inhibitor signal
    diffusion_pixels = um_to_pixels(diffusion_um, image_width_px, film_format_mm)
    if diffusion_pixels > EPSILON:
        # We blur the density signal, which is proportional to the amount of inhibitor released
        inhibitor_signal_diffused = gaussian_filter(naive_density_rgb, sigma=(diffusion_pixels, diffusion_pixels, 0))
    else:
        inhibitor_signal_diffused = naive_density_rgb

    # 2. Apply inter-layer crosstalk
    # inhibitor_signal has shape (H, W, 3), inhibitor_matrix has shape (3, 3)
    # einsum calculates the total inhibition for each destination layer from all source layers
    inhibitor_effect = np.einsum('...s, sm -> ...m', inhibitor_signal_diffused, inhibitor_matrix)
    
    # 3. Modify the original log exposure
    modified_log_exposure = log_exposure_rgb - inhibitor_effect
    return modified_log_exposure


def apply_hd_curves(
    log_exposure_rgb,
    processing: Processing,
    hd_curve: List[HDCurvePoint],
    middle_gray_logE: float,
) -> np.ndarray:
    """Applies H&D curves to log exposure values."""
    density_rgb = np.zeros_like(log_exposure_rgb)
    gamma_factors = [
        processing.gamma.r_factor,
        processing.gamma.g_factor,
        processing.gamma.b_factor,
    ]
    balance_shifts = [
        processing.balance.r_shift,
        processing.balance.g_shift,
        processing.balance.b_shift,
    ]

    min_logE = hd_curve[0].d
    max_logE = hd_curve[-1].d
    min_densities = [hd_curve[0].g, hd_curve[0].g, hd_curve[0].g]
    max_densities = [hd_curve[-1].g, hd_curve[-1].g, hd_curve[-1].g]

    for i, channel in enumerate(["R", "G", "B"]):
        # Apply gamma factor (affects contrast by scaling log exposure input)
        # Handle potential division by zero if gamma factor is 0
        gamma_factor = gamma_factors[i]
        if abs(gamma_factor) < EPSILON:
            print(
                f"Warning: Gamma factor for channel {channel} is near zero. Clamping to {EPSILON}.",
                file=sys.stderr,
            )
            gamma_factor = EPSILON if gamma_factor >= 0 else -EPSILON

        # Adjust log exposure relative to middle gray before applying gamma
        log_exposure_adjusted = (
            middle_gray_logE
            + (log_exposure_rgb[..., i] - middle_gray_logE) / gamma_factor
        )

        # Clamp adjusted exposure to the range defined in the H&D data before interpolation
        log_exposure_clamped = np.clip(log_exposure_adjusted, min_logE, max_logE)

        # Create interpolation function for the current channel
        if channel == "R":
            interp_func = interp1d(
                [d.d for d in hd_curve],
                [d.r for d in hd_curve],
                kind="linear",  # Linear interpolation is common, could be 'cubic'
                bounds_error=False,  # Allows extrapolation, but we clamp manually below
                fill_value=(min_densities[i], max_densities[i]),  # type: ignore
            )
        elif channel == "G":
            interp_func = interp1d(
                [d.d for d in hd_curve],
                [d.g for d in hd_curve],
                kind="linear",  # Linear interpolation is common, could be 'cubic'
                bounds_error=False,  # Allows extrapolation, but we clamp manually below
                fill_value=(min_densities[i], max_densities[i]),  # type: ignore
            )
        else:
            interp_func = interp1d(
                [d.d for d in hd_curve],
                [d.b for d in hd_curve],
                kind="linear",  # Linear interpolation is common, could be 'cubic'
                bounds_error=False,  # Allows extrapolation, but we clamp manually below
                fill_value=(min_densities[i], max_densities[i]),  # type: ignore
            )

        # Apply interpolation
        density = interp_func(log_exposure_clamped)

        # Apply density balance shift (additive density offset)
        density += balance_shifts[i]

        density_rgb[..., i] = np.maximum(density, 0)  # Ensure density is non-negative

    return density_rgb


def apply_saturation_rgb(image_linear, saturation_factor):
    """Adjusts saturation directly in RGB space."""
    if saturation_factor == 1.0:
        return image_linear

    # Luminance weights for sRGB primaries (Rec.709)
    luminance = (
        0.2126 * image_linear[..., 0]
        + 0.7152 * image_linear[..., 1]
        + 0.0722 * image_linear[..., 2]
    )

    # Expand luminance to 3 channels for broadcasting
    luminance_rgb = np.expand_dims(luminance, axis=-1)

    # Apply saturation: Lerp between luminance (grayscale) and original color
    saturated_image = luminance_rgb + saturation_factor * (image_linear - luminance_rgb)

    # Clip results to valid range (important after saturation boost)
    return np.clip(saturated_image, 0.0, 1.0)


def apply_spatial_effects(
    image,
    film_format_mm,
    interlayerData: Interlayer,
    halationData: Halation,
    image_width_px,
):
    """Applies diffusion blur and halation."""
    # Combine diffusion effects (assuming they add quadratically in terms of sigma)
    softening_diffusion_um = interlayerData.diffusion_um
    if softening_diffusion_um > EPSILON:
        sigma_pixels_diffusion = um_to_pixels(softening_diffusion_um, image_width_px, film_format_mm)
        if sigma_pixels_diffusion > EPSILON:
            print(f"Applying interlayer diffusion blur: sigma={sigma_pixels_diffusion:.2f} pixels ({softening_diffusion_um:.1f} um)")
            image = gaussian_filter(image, sigma=[sigma_pixels_diffusion, sigma_pixels_diffusion, 0], mode="nearest")
            image = np.clip(image, 0.0, 1.0)

    # --- 2. Apply Halation ---
    # This simulates light scattering back through the emulsion
    halation_applied = True
    blurred_image_halation = np.copy(
        image
    )  # Start with potentially diffusion-blurred image

    strengths = [
        halationData.strength.r,
        halationData.strength.g,
        halationData.strength.b,
    ]
    sizes_um = [
        halationData.size_um.r,
        halationData.size_um.g,
        halationData.size_um.b,
    ]

    for i in range(3):
        strength = strengths[i]
        size_um = sizes_um[i]
        if strength > EPSILON and size_um > EPSILON:
            sigma_pixels_halation = um_to_pixels(
                size_um, image_width_px, film_format_mm
            )
            if sigma_pixels_halation > EPSILON:
                halation_applied = True
                print(
                    f"Applying halation blur (Channel {i}): sigma={sigma_pixels_halation:.2f} pixels ({size_um:.1f} um), strength={strength:.3f}"
                )
                # Blur only the current channel for halation effect
                channel_blurred = gaussian_filter(
                    image[..., i], sigma=sigma_pixels_halation, mode="nearest"
                )
                # Add the blurred channel back, weighted by strength
                blurred_image_halation[..., i] = (
                    image[..., i] * (1.0 - strength) + channel_blurred * strength
                )

    if halation_applied:
        # Clip final result after halation
        image = np.clip(blurred_image_halation, 0.0, 1.0)

    return image


def calculate_film_log_exposure(image_linear_srgb, # Assuming input is converted to linear sRGB
                                spectral_data: list[SpectralSensitivityCurvePoint],
                                ref_illuminant_spd, # e.g., colour.SDS_ILLUMINANTS['D65']
                                middle_gray_logE,
                                EPSILON=1e-10):
    """
    Calculates the effective log exposure for each film layer based on spectral sensitivities.
    This version correctly maps film layers (RGB) to dye sensitivities (CMY)
    and calibrates exposure without distorting color balance.
    """
    # --- 1. Define a common spectral shape and align all data ---
    common_shape = colour.SpectralShape(380, 780, 5)
    common_wavelengths = common_shape.wavelengths
    
    # --- THE FIX IS HERE: Correct mapping of film layers to sensitivities ---
    # The 'R' layer of our model corresponds to the Cyan dye ('c') sensitivity curve.
    # The 'G' layer of our model corresponds to the Magenta dye ('m') sensitivity curve.
    # The 'B' layer of our model corresponds to the Yellow dye ('y') sensitivity curve.
    film_sens_R = interp1d([p.wavelength for p in spectral_data], [p.c for p in spectral_data], bounds_error=False, fill_value=0)(common_wavelengths)
    film_sens_G = interp1d([p.wavelength for p in spectral_data], [p.m for p in spectral_data], bounds_error=False, fill_value=0)(common_wavelengths)
    film_sens_B = interp1d([p.wavelength for p in spectral_data], [p.y for p in spectral_data], bounds_error=False, fill_value=0)(common_wavelengths)
    
    film_spectral_sensitivities = np.stack([film_sens_R, film_sens_G, film_sens_B], axis=-1)
    print(f"Film spectral sensitivities shape: {film_spectral_sensitivities.shape}")

    # Align Reference Illuminant to our common shape
    illuminant_aligned = ref_illuminant_spd.copy().align(common_shape)

    # --- 2. Use Mallett (2019) to get spectral reflectance from sRGB input ---
    # Get the three sRGB spectral primary basis functions
    mallett_basis_sds = colour.recovery.MSDS_BASIS_FUNCTIONS_sRGB_MALLETT2019
    mallett_basis_aligned = mallett_basis_sds.copy().align(common_shape)
    print(f"Mallett basis shape: {mallett_basis_aligned.shape}")

    # This is the core of Mallett's method: the reflectance spectrum of any sRGB color
    # is a linear combination of these three basis spectra, weighted by the linear sRGB values.
    # S(lambda) = R * Sr(lambda) + G * Sg(lambda) + B * Sb(lambda)
    spectral_reflectance = np.einsum(
        '...c, kc -> ...k',  # c: sRGB channels, k: wavelengths
        image_linear_srgb, mallett_basis_aligned.values
    )
    print(f"Spectral reflectance shape: {spectral_reflectance.shape}")
    
    # --- 3. Calculate Scene Light and Film Exposure ---
    # The light hitting the film is the spectral reflectance multiplied by the scene illuminant.
    light_from_scene = spectral_reflectance * illuminant_aligned.values
    print(f"Light from scene shape: {light_from_scene.shape}")
    
    # The exposure on each film layer is the integral of scene light * layer sensitivity.
    # The einsum here performs that multiplication and summation (integration) in one step.
    film_exposure_values = np.einsum(
        '...k, ks -> ...s', # k: wavelengths, s: film layers
        light_from_scene, film_spectral_sensitivities
    )
    print(f"Film exposure values shape: {film_exposure_values.shape}")

    # --- 4. Calibrate Exposure Level (The Right Way) ---
    # We now need to anchor this exposure to the datasheet's middle_gray_logE.
    # First, find the exposure produced by a perfect 18.4% gray card.
    gray_srgb_linear = np.array([0.184, 0.184, 0.184])
    gray_reflectance = np.einsum('c, kc -> k', gray_srgb_linear, mallett_basis_aligned.values)
    gray_light = gray_reflectance * illuminant_aligned.values
    exposure_18_gray_film = np.einsum('k, ks -> s', gray_light, film_spectral_sensitivities)
    print(f"Exposure for 18.4% gray card shape: {exposure_18_gray_film.shape}")

    # The datasheet says that a middle gray exposure should result in a log value of -1.44 (for Portra).
    # Our "green" channel is the usual reference for overall brightness.
    # Let's find out what log exposure our 18.4% gray card *actually* produced on the green layer.
    log_exposure_of_gray_on_green_layer = np.log10(exposure_18_gray_film[1] + EPSILON)
    
    # Now, calculate the difference (the "shift") needed to match the datasheet.
    log_shift = middle_gray_logE - log_exposure_of_gray_on_green_layer
    print(f"Log shift to match middle gray: {log_shift:.3f}")

    # --- 5. Final Conversion to Log Exposure ---
    # Apply this single brightness shift to all three channels. This preserves the
    # film's inherent color balance while correctly setting the exposure level.
    log_exposure_rgb = np.log10(film_exposure_values + EPSILON) + log_shift
    print(f"Log exposure RGB shape: {log_exposure_rgb.shape}")

    return log_exposure_rgb


def main():
    parser = argparse.ArgumentParser(
        description="Simulate film stock color characteristics using a datasheet JSON."
    )
    parser.add_argument(
        "input_image",
        help="Path to the input RGB image (e.g., PNG, TIFF). Assumed linear RGB.",
    )
    parser.add_argument("datasheet_json", help="Path to the film datasheet JSON file.")
    parser.add_argument("output_image", help="Path to save the output emulated image.")

    args = parser.parse_args()

    # --- Load Datasheet ---
    print(f"Loading datasheet: {args.datasheet_json}")
    datasheet: FilmDatasheet | None = parse_datasheet_json(args.datasheet_json)
    if datasheet is None:
        sys.exit(1)
    print(
        f"Simulating: {datasheet.info.name} ({datasheet.info.format_mm}mm) (v{datasheet.info.version})\n\t{datasheet.info.description}"
    )

    print("Pre-calculating DIR coupler effects...")
    inhibitor_matrix = compute_inhibitor_matrix(datasheet.properties.couplers.dir_amount_rgb, datasheet.properties.couplers.dir_diffusion_interlayer)
    print("Inhibitor Matrix:\n", inhibitor_matrix)
    uncoupled_hd_curve = compute_uncoupled_hd_curves(datasheet.properties.curves.hd, inhibitor_matrix)
    print(f"Successfully computed {len(uncoupled_hd_curve)} points for the uncoupled H&D curve.")

    import pprint
    pprint.pp(datasheet)

    # --- Load Input Image ---
    print(f"Loading input image: {args.input_image}")
    try:
        # For DNG files, force reading the raw image data, not the embedded thumbnail
        if (
            args.input_image.lower().endswith(".dng")
            or args.input_image.lower().endswith(".raw")
            or args.input_image.lower().endswith(".arw")
        ):
            print("Detected Camera RAW file, reading raw image data...")
            with rawpy.imread(args.input_image) as raw:
                image_raw = raw.postprocess(
                    demosaic_algorithm=rawpy.DemosaicAlgorithm.AHD,  # type: ignore
                    output_bps=16,  # Use 16-bit output for better precision
                    use_camera_wb=True,  # Use camera white balance
                    no_auto_bright=True,  # Disable auto brightness adjustment
                    output_color=rawpy.ColorSpace.sRGB,  # type: ignore
                )
            # If the image has more than 3 channels, try to select the first 3 (RGB)
            if image_raw.ndim == 3 and image_raw.shape[-1] > 3:
                image_raw = image_raw[..., :3]
        else:
            image_raw = iio.imread(args.input_image)
    except FileNotFoundError:
        print(
            f"Error: Input image file not found at {args.input_image}", file=sys.stderr
        )
        sys.exit(1)
    except Exception as e:
        print(f"Error reading input image: {e}", file=sys.stderr)
        sys.exit(1)

    # Check if image is likely sRGB (8-bit or 16-bit integer types are usually sRGB)
    is_srgb = image_raw.dtype in (np.uint8, np.uint16)
    if is_srgb:
        print("Input image appears to be sRGB. Linearizing...")
        # Convert to float in [0,1]
        if image_raw.dtype == np.uint8:
            image_float = image_raw.astype(np.float64) / 255.0
        elif image_raw.dtype == np.uint16:
            image_float = image_raw.astype(np.float64) / 65535.0
        else:
            image_float = image_raw.astype(np.float64)

        # sRGB to linear conversion
        def srgb_to_linear(c):
            c = np.clip(c, 0.0, 1.0)
            return np.where(c <= 0.04045, c / 12.92, ((c + 0.055) / 1.055) ** 2.4)

        image_raw = srgb_to_linear(image_float)
    else:
        print("Input image is assumed to be linear RGB.")

    # --- Prepare Image Data ---
    # Convert to float64 for precision, handle different input types
    if image_raw.dtype == np.uint8:
        image_linear = image_raw.astype(np.float64) / 255.0
    elif image_raw.dtype == np.uint16:
        image_linear = image_raw.astype(np.float64) / 65535.0
    elif image_raw.dtype == np.float32:
        image_linear = image_raw.astype(np.float64)
    elif image_raw.dtype == np.float64:
        image_linear = image_raw
    else:
        print(f"Error: Unsupported image data type: {image_raw.dtype}", file=sys.stderr)
        sys.exit(1)

    # Discard alpha channel if present
    if image_linear.shape[-1] == 4:
        print("Discarding alpha channel.")
        image_linear = image_linear[..., :3]
    elif image_linear.shape[-1] != 3:
        print(
            f"Error: Input image must be RGB (shape {image_linear.shape} not supported).",
            file=sys.stderr,
        )
        sys.exit(1)

    # Ensure input is non-negative
    image_linear = np.maximum(image_linear, 0.0)
    print(f"Input image dimensions: {image_linear.shape}")
    image_width_px = image_linear.shape[1]

    # --- Pipeline Steps ---

    # 1. Convert Linear RGB to Log Exposure (LogE)
    # Map linear 0.18 to the specified middle_gray_logE from the datasheet
    print("Converting linear RGB to Log Exposure...")
    middle_gray_logE = float(datasheet.properties.calibration.middle_gray_logE)
    # Add epsilon inside log10 to handle pure black pixels
    log_exposure_rgb = calculate_film_log_exposure(image_linear, 
                                                   datasheet.properties.curves.spectral_sensitivity,
                                                   SDS_ILLUMINANTS['D65'],  # Use D65 as reference illuminant
                                                   middle_gray_logE, EPSILON)
    

    print("Applying DIR coupler simulation...")
    naive_density_rgb = apply_hd_curves(log_exposure_rgb, datasheet.processing, datasheet.properties.curves.hd, middle_gray_logE)
    # 2b. Use this density to calculate the inhibitor effect and modify the log exposure
    modified_log_exposure_rgb = apply_dir_coupler_simulation(
        log_exposure_rgb,
        naive_density_rgb,
        inhibitor_matrix,
        datasheet.properties.couplers.dir_diffusion_um,
        datasheet.info.format_mm,
        image_width_px
    )


    # 2. Apply H&D Curves (Tonal Mapping + Balance Shifts + Gamma/Contrast)
    print("Applying H&D curves...")
    density_rgb = apply_hd_curves(
        modified_log_exposure_rgb,
        datasheet.processing,
        datasheet.properties.curves.hd,
        middle_gray_logE,
    )

    # 3. Convert Density back to Linear Transmittance
    # Higher density means lower transmittance
    print("Converting density to linear transmittance...")
    # Add density epsilon? Usually density floor (Dmin) handles this.
    linear_transmittance = 10.0 ** (-density_rgb)
    # Normalize transmittance? Optional, assumes Dmin corresponds roughly to max transmittance 1.0
    # Could normalize relative to Dmin from the curves if needed.
    # linear_transmittance = linear_transmittance / (10.0**(-np.array([hd_data['Density_R'][0], hd_data['Density_G'][0], hd_data['Density_B'][0]])))
    linear_transmittance = np.clip(linear_transmittance, 0.0, 1.0)

    # 4. Apply Spatial Effects (Diffusion Blur, Halation)
    print("Applying spatial effects (diffusion, halation)...")
    # Apply these effects in the linear domain
    linear_post_spatial = apply_spatial_effects(
        linear_transmittance,
        datasheet.info.format_mm,
        datasheet.properties.interlayer,
        datasheet.properties.halation,
        image_width_px,
    )

    # 5. Apply Saturation Adjustment (Approximating Coupler Effects)
    # Assuming coupler_amount directly scales saturation factor.
    # Values > 1 increase saturation, < 1 decrease.
    linear_post_saturation = apply_saturation_rgb(linear_post_spatial, datasheet.properties.couplers.saturation_amount)

    # --- Final Output Conversion ---
    print("Converting to output format...")
    # Clip final result and convert to uint8
    if args.output_image.lower().endswith(".tiff"):
        output_image_uint8 = (
            np.clip(linear_post_saturation, 0.0, 1.0) * 65535.0
        ).astype(np.uint16)
    else:
        output_image_uint8 = (np.clip(linear_post_saturation, 0.0, 1.0) * 255.0).astype(
            np.uint8
        )

    # --- Save Output Image ---
    print(f"Saving output image: {args.output_image}")
    try:
        if args.output_image.lower().endswith((".tiff", ".tif")):
            # Use imageio for standard formats
            iio.imwrite(args.output_image, output_image_uint8)
        elif args.output_image.lower().endswith(".png"):
            iio.imwrite(args.output_image, output_image_uint8, format="PNG")
        else:
            iio.imwrite(args.output_image, output_image_uint8, quality=95)
        print("Done.")
    except Exception as e:
        print(f"Error writing output image: {e}", file=sys.stderr)
        sys.exit(1)


if __name__ == "__main__":
    main()