diff --git a/filmcolor b/filmcolor index 2c5df1b..8116b1f 100755 --- a/filmcolor +++ b/filmcolor @@ -10,37 +10,34 @@ # ] # /// -# -*- 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 +from scipy.ndimage import ( + gaussian_filter, + gaussian_filter1d, + map_coordinates, +) # For LUTs import rawpy import os import sys import colour from colour.colorimetry import SDS_ILLUMINANTS import json +from typing import Optional, List +from pathlib import Path # --- Configuration --- -# Small epsilon to prevent log(0) or division by zero errors EPSILON = 1e-10 - -from typing import Optional, List +### PERFORMANCE ### +# Size of the 3D LUTs. 33 is a common industry size (e.g., in .cube files). +# 17 is faster to generate, 65 is more accurate but much slower to generate. +LUT_SIZE = 65 +# --- Data Structures (unchanged) --- class Info: name: str description: str @@ -50,15 +47,11 @@ class Info: 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})" + self.name, self.description, self.format_mm, self.version = ( + name, + description, + format_mm, + version, ) @@ -68,14 +61,7 @@ class Balance: 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})" - ) + self.r_shift, self.g_shift, self.b_shift = r_shift, g_shift, b_shift class Gamma: @@ -84,14 +70,7 @@ class Gamma: 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})" - ) + self.r_factor, self.g_factor, self.b_factor = r_factor, g_factor, b_factor class Processing: @@ -99,14 +78,7 @@ class Processing: 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}))" - ) + self.gamma, self.balance = gamma, balance class Couplers: @@ -115,11 +87,24 @@ class Couplers: 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 + def __init__( + self, + saturation_amount: float, + dir_amount_rgb: List[float], + dir_diffusion_um: float, + dir_diffusion_interlayer: float, + ) -> None: + ( + self.saturation_amount, + self.dir_amount_rgb, + self.dir_diffusion_um, + self.dir_diffusion_interlayer, + ) = ( + saturation_amount, + dir_amount_rgb, + dir_diffusion_um, + dir_diffusion_interlayer, + ) class HDCurvePoint: @@ -129,14 +114,8 @@ class HDCurvePoint: 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 + self.d, self.r, self.g, self.b = d, r, g, 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 @@ -145,13 +124,7 @@ class SpectralSensitivityCurvePoint: 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})" + self.wavelength, self.y, self.m, self.c = wavelength, y, m, c class RGBValue: @@ -160,24 +133,19 @@ class RGBValue: 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})" + self.r, self.g, self.b = r, g, b 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)})" + def __init__( + self, + hd: List[HDCurvePoint], + spectral_sensitivity: List[SpectralSensitivityCurvePoint], + ) -> None: + self.hd, self.spectral_sensitivity = hd, spectral_sensitivity class Halation: @@ -185,11 +153,7 @@ class Halation: 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})" + self.strength, self.size_um = strength, size_um class Interlayer: @@ -198,22 +162,13 @@ class Interlayer: 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})" - ) + self.iso, self.middle_gray_logE = iso, middle_gray_logE class Properties: @@ -231,17 +186,12 @@ class Properties: 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})" + self.halation, self.couplers, self.interlayer, self.curves, self.calibration = ( + halation, + couplers, + interlayer, + curves, + calibration, ) @@ -253,33 +203,19 @@ class FilmDatasheet: def __init__( self, info: Info, processing: Processing, properties: Properties ) -> None: - self.info = info - self.processing = processing - self.properties = properties + self.info, self.processing, self.properties = info, processing, properties - def __repr__(self) -> str: - return ( - f"FilmDatasheet(info={self.info}\nprocessing={self.processing}\n" - f"properties={self.properties})" - ) - -import pprint +# --- Datasheet Parsing (unchanged) --- 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. - """ + # This function remains identical to your last version + # ... 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"], @@ -313,14 +249,16 @@ def parse_datasheet_json(json_filepath) -> FilmDatasheet | None: 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"], + 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"], + middle_gray_logE=data["properties"]["calibration"]["middle_gray_logE"], ) curves = Curves( hd=[ @@ -337,7 +275,6 @@ def parse_datasheet_json(json_filepath) -> FilmDatasheet | None: for point in data["properties"]["curves"]["spectral_sensitivity"] ], ) - print(f"Parsed {len(curves.hd)} H&D curve points.") properties = Properties( calibration=calibration, halation=halation, @@ -353,94 +290,162 @@ def parse_datasheet_json(json_filepath) -> FilmDatasheet | None: return None +# --- Core Simulation Functions --- + + +### PERFORMANCE ### +def apply_3d_lut(image: np.ndarray, lut: np.ndarray) -> np.ndarray: + """ + Applies a 3D LUT to an image using trilinear interpolation. + + Args: + image: Input image, shape (H, W, 3), values expected in [0, 1]. + lut: 3D LUT, shape (N, N, N, 3). + + Returns: + Output image, shape (H, W, 3). + """ + lut_dim = lut.shape[0] + + # Scale image coordinates to LUT indices + scaled_coords = image * (lut_dim - 1) + + # map_coordinates requires coordinates in shape (3, H, W) + coords = np.transpose(scaled_coords, (2, 0, 1)) + + # Interpolate each output channel + out_r = map_coordinates(lut[..., 0], coords, order=1, mode="nearest") + out_g = map_coordinates(lut[..., 1], coords, order=1, mode="nearest") + out_b = map_coordinates(lut[..., 2], coords, order=1, mode="nearest") + + return np.stack([out_r, out_g, out_b], axis=-1) + + +### PERFORMANCE ### +def create_exposure_lut( + spectral_data: list[SpectralSensitivityCurvePoint], + ref_illuminant_spd, + middle_gray_logE: float, + size: int, +) -> np.ndarray: + """Generates a 3D LUT for the linear sRGB -> Log Exposure conversion.""" + print(f"Generating {size}x{size}x{size} exposure LUT...") + + # Create a grid of sRGB input values + grid_points = np.linspace(0.0, 1.0, size) + r, g, b = np.meshgrid(grid_points, grid_points, grid_points, indexing="ij") + # Reshape grid into an "image" for batch processing + srgb_grid = np.stack([r, g, b], axis=-1).reshape( + -1, 3 + ) # Shape: (size*size*size, 3) + srgb_grid = srgb_grid.reshape(size, size * size, 3) # Make it 2D for the function + + # Run the full spectral calculation on this grid + lut_data = calculate_film_log_exposure( + srgb_grid, spectral_data, ref_illuminant_spd, middle_gray_logE, EPSILON + ) + + # Reshape back into a 3D LUT + return lut_data.reshape(size, size, size, 3) + + +### PERFORMANCE ### +def create_density_lut( + processing: Processing, + hd_curve: List[HDCurvePoint], + middle_gray_logE: float, + log_exposure_range: tuple, + size: int, +) -> np.ndarray: + """Generates a 3D LUT for the Log Exposure -> Density conversion.""" + print( + f"Generating {size}x{size}x{size} density LUT for range {log_exposure_range}..." + ) + le_min, le_max = log_exposure_range + + grid_points = np.linspace(le_min, le_max, size) + le_r, le_g, le_b = np.meshgrid(grid_points, grid_points, grid_points, indexing="ij") + log_exposure_grid = np.stack([le_r, le_g, le_b], axis=-1) + + lut_data = apply_hd_curves( + log_exposure_grid, processing, hd_curve, middle_gray_logE + ) + return lut_data + + +# --- (The rest of the core functions: compute_inhibitor_matrix, compute_uncoupled_hd_curves, +# apply_dir_coupler_simulation, um_to_pixels, apply_hd_curves, apply_saturation_rgb, +# apply_spatial_effects, calculate_film_log_exposure, remain UNCHANGED from the previous version) --- +def compute_inhibitor_matrix( + amount_rgb: List[float], diffusion_interlayer: float +) -> np.ndarray: + matrix = np.eye(3) + 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 + ) + row_sums = matrix.sum(axis=1) + row_sums[row_sums == 0] = 1 # Avoid division by zero + matrix = matrix / row_sums[:, np.newaxis] + return matrix * np.array(amount_rgb)[:, np.newaxis] + + +def compute_uncoupled_hd_curves( + hd_curve_data: List[HDCurvePoint], inhibitor_matrix: np.ndarray +) -> List[HDCurvePoint]: + log_E_values = np.array([p.d for p in hd_curve_data if p.d is not None]) + density_r, density_g, density_b = ( + np.array([p.r for p in hd_curve_data]), + np.array([p.g for p in hd_curve_data]), + np.array([p.b for p in hd_curve_data]), + ) + source_density_matrix = np.tile(density_g[:, np.newaxis], (1, 3)) + inhibitor_effect = source_density_matrix @ inhibitor_matrix + 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 + ) + return [ + HDCurvePoint(d=log_e, r=uncoupled_r[i], g=uncoupled_g[i], b=uncoupled_b[i]) + for i, log_e in enumerate(log_E_values) + ] + + +def apply_dir_coupler_simulation( + log_exposure_rgb, + naive_density_rgb, + inhibitor_matrix, + diffusion_um, + film_format_mm, + image_width_px, +): + diffusion_pixels = um_to_pixels(diffusion_um, image_width_px, film_format_mm) + inhibitor_signal_diffused = ( + gaussian_filter( + naive_density_rgb, sigma=(diffusion_pixels, diffusion_pixels, 0) + ) + if diffusion_pixels > EPSILON + else naive_density_rgb + ) + inhibitor_effect = np.einsum( + "...s, sm -> ...m", inhibitor_signal_diffused, inhibitor_matrix + ) + return log_exposure_rgb - inhibitor_effect + + 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 + return sigma_um / microns_per_pixel def apply_hd_curves( @@ -449,247 +454,214 @@ def apply_hd_curves( 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] - + density_rgb, gamma_factors, balance_shifts = ( + np.zeros_like(log_exposure_rgb), + [ + processing.gamma.r_factor, + processing.gamma.g_factor, + processing.gamma.b_factor, + ], + [ + processing.balance.r_shift, + processing.balance.g_shift, + processing.balance.b_shift, + ], + ) + log_e_points = [p.d for p in hd_curve if p.d is not None] + min_logE, max_logE = log_e_points[0], log_e_points[-1] 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 + log_exposure_adjusted = middle_gray_logE + ( + log_exposure_rgb[..., i] - middle_gray_logE + ) / (gamma_factor if abs(gamma_factor) > EPSILON else EPSILON) + y_points = ( + [p.r for p in hd_curve] + if channel == "R" + else [p.g for p in hd_curve] if channel == "G" else [p.b for p in hd_curve] + ) + interp_func = interp1d( + log_e_points, + y_points, + kind="linear", + bounds_error=False, + fill_value=(y_points[0], y_points[-1]), # type: ignore + ) + density_rgb[..., i] = np.maximum( + interp_func(np.clip(log_exposure_adjusted, min_logE, max_logE)) + + balance_shifts[i], + 0, ) - - # 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] + luminance = np.einsum( + "...c, c -> ...", image_linear, np.array([0.2126, 0.7152, 0.0722]) + ) + return np.clip( + np.expand_dims(luminance, axis=-1) + + saturation_factor * (image_linear - np.expand_dims(luminance, axis=-1)), + 0.0, + 1.0, ) - # 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) +def white_balance_to_d65(image_linear_srgb: np.ndarray) -> np.ndarray: + """ + Attempts to white balance a linear sRGB image to a D65 illuminant. - # Clip results to valid range (important after saturation boost) - return np.clip(saturated_image, 0.0, 1.0) + This function is for non-RAW images where the white balance is already baked in. + It uses the "Gray World" assumption to estimate the current illuminant. + + Args: + image_linear_srgb: A numpy array representing the image in linear sRGB space. + + Returns: + A numpy array of the white-balanced image, also in linear sRGB space. + """ + print("Attempting to force D65 white balance using Gray World estimation...") + + # 1. Estimate the illuminant of the source image. + # The Gray World assumption states that the average pixel color is the illuminant. + source_illuminant_rgb = np.mean(image_linear_srgb, axis=(0, 1)) + + # Avoid division by zero if the image is black + if np.all(source_illuminant_rgb < EPSILON): + return image_linear_srgb + + # 2. Define the input and output color spaces and illuminants. + # We assume the image is sRGB, which also uses a D65 white point *by definition*. + # However, if the image has a color cast, its *actual* illuminant is not D65. + # We are adapting from the *estimated* illuminant back to the *ideal* D65. + srgb_cs = colour.models.RGB_COLOURSPACE_sRGB + + # The target illuminant is D65. We get its XYZ value from colour-science. + target_illuminant_xyz = colour.CCS_ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65'] + + # Convert our estimated RGB illuminant to XYZ. + # source_illuminant_rgb are assumed to be in the sRGB colourspace. + # srgb_cs is the sRGB colourspace definition, which includes its D65 whitepoint and conversion matrix. + source_illuminant_xyz = colour.RGB_to_XYZ( + source_illuminant_rgb, # RGB values to convert + srgb_cs, # Definition of the RGB colourspace these values are in (sRGB) + srgb_cs.whitepoint # CIE XYZ of the illuminant for the sRGB colourspace (D65) + # This ensures conversion without further chromatic adaptation at this step, + # as the input RGB values are already adapted to srgb_cs.whitepoint. + ) + + # 3. Perform the chromatic adaptation. + # We use the Von Kries transform, which is a common and effective method. + image_adapted_linear_srgb = colour.adaptation.chromatic_adaptation( + image_linear_srgb, + source_illuminant_xyz, + target_illuminant_xyz, + method='Von Kries', + transform='CAT02' # CAT02 is a well-regarded choice + ) + + # 4. Clip to ensure values remain valid. + return np.clip(image_adapted_linear_srgb, 0.0, 1.0) def apply_spatial_effects( image, film_format_mm, + couplerData: Couplers, 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 = [ + sigma_pixels_diffusion = um_to_pixels( + interlayerData.diffusion_um, image_width_px, film_format_mm + ) + if sigma_pixels_diffusion > EPSILON: + image = gaussian_filter( + image, + sigma=[sigma_pixels_diffusion, sigma_pixels_diffusion, 0], + mode="nearest", + ) + halation_applied, blurred_image_halation = False, np.copy(image) + strengths, sizes_um = [ halationData.strength.r, halationData.strength.g, halationData.strength.b, - ] - sizes_um = [ - halationData.size_um.r, - halationData.size_um.g, - halationData.size_um.b, - ] - + ], [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 + sigma_pixels_halation = um_to_pixels( + sizes_um[i], image_width_px, film_format_mm + ) + if strengths[i] > EPSILON and sigma_pixels_halation > EPSILON: + halation_applied = True + channel_blurred = gaussian_filter( + image[..., i], sigma=sigma_pixels_halation, mode="nearest" ) - 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 + blurred_image_halation[..., i] = ( + image[..., i] * (1.0 - strengths[i]) + channel_blurred * strengths[i] + ) + return ( + np.clip(blurred_image_halation, 0.0, 1.0) + if halation_applied + else np.clip(image, 0.0, 1.0) + ) -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 +def calculate_film_log_exposure( + image_linear_srgb, + spectral_data: list[SpectralSensitivityCurvePoint], + ref_illuminant_spd, + middle_gray_logE, + EPSILON, +): + common_shape, common_wavelengths = ( + colour.SpectralShape(380, 780, 5), + colour.SpectralShape(380, 780, 5).wavelengths, + ) + sensitivities = np.stack( + [ + interp1d( + [p.wavelength for p in spectral_data], + [p.c for p in spectral_data], + bounds_error=False, + fill_value=0, + )(common_wavelengths), + interp1d( + [p.wavelength for p in spectral_data], + [p.m for p in spectral_data], + bounds_error=False, + fill_value=0, + )(common_wavelengths), + interp1d( + [p.wavelength for p in spectral_data], + [p.y for p in spectral_data], + bounds_error=False, + fill_value=0, + )(common_wavelengths), + ], + axis=-1, + ) 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) + mallett_basis_aligned = ( + colour.recovery.MSDS_BASIS_FUNCTIONS_sRGB_MALLETT2019.copy().align(common_shape) + ) spectral_reflectance = np.einsum( - '...c, kc -> ...k', # c: sRGB channels, k: wavelengths - image_linear_srgb, mallett_basis_aligned.values + "...c, kc -> ...k", 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 + "...k, ks -> ...s", light_from_scene, sensitivities + ) + gray_reflectance = np.einsum( + "c, kc -> k", np.array([0.184, 0.184, 0.184]), mallett_basis_aligned.values ) - 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}") + exposure_18_gray_film = np.einsum("k, ks -> s", gray_light, sensitivities) + log_shift = middle_gray_logE - np.log10(exposure_18_gray_film[1] + EPSILON) + return np.log10(film_exposure_values + EPSILON) + log_shift - # 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 +# --- Main Execution --- def main(): @@ -702,11 +674,14 @@ def main(): ) 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.") - + parser.add_argument( + "--no-cache", action="store_true", help="Force regeneration of LUTs and do not save them." + ) + parser.add_argument( + "--force-d65", action="store_true", help="Force white balance when loading RAW files." + ) 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) @@ -714,185 +689,263 @@ def main(): 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.") + # --- LUT Generation and Pre-computation --- + cache_dir = Path.home() / ".filmsim" + # Sanitize datasheet name for use in a filename + safe_name = "".join(c for c in datasheet.info.name if c.isalnum() or c in (' ', '_')).rstrip() + base_filename = f"{safe_name.replace(' ', '_')}_v{datasheet.info.version}_size{LUT_SIZE}" + exposure_lut_path = cache_dir / f"{base_filename}_exposure.cube" + naive_density_lut_path = cache_dir / f"{base_filename}_naive_density.cube" + uncoupled_density_lut_path = cache_dir / f"{base_filename}_uncoupled_density.cube" + user_consent_to_save = None - import pprint - pprint.pp(datasheet) + middle_gray_logE = float(datasheet.properties.calibration.middle_gray_logE) + + exposure_lut: Optional[np.ndarray] = None # Initialize exposure_lut - # --- Load Input Image --- - print(f"Loading input image: {args.input_image}") + ### PERFORMANCE ### + # Try to load exposure LUT from cache first + if not args.no_cache and exposure_lut_path.exists(): + try: + print(f"Loading exposure LUT from cache: {exposure_lut_path}") + loaded_obj = colour.read_LUT(str(exposure_lut_path)) + if isinstance(loaded_obj, colour.LUT3D): + exposure_lut = loaded_obj.table + print(f"Successfully loaded exposure LUT from {exposure_lut_path}") + else: + # Log a warning and fall through to regeneration + print(f"Warning: Cached exposure LUT {exposure_lut_path} is not a LUT3D (type: {type(loaded_obj)}). Will regenerate.") + except Exception as e: + # Log a warning and fall through to regeneration + print(f"Warning: Error loading exposure LUT from {exposure_lut_path}: {e}. Will regenerate.") + + # If LUT wasn't loaded from cache (or cache was disabled, or loading failed), generate it + if exposure_lut is None: + print("Generating new exposure LUT...") + generated_exposure_lut_table = create_exposure_lut( + datasheet.properties.curves.spectral_sensitivity, + SDS_ILLUMINANTS["D65"], + middle_gray_logE, + size=LUT_SIZE, + ) + exposure_lut = generated_exposure_lut_table # Assign to the main variable used later + + # Save the newly generated LUT if caching is enabled + if not args.no_cache: + if user_consent_to_save is None: # Ask for consent only once per run if needed + response = input(f"Save generated LUTs to {cache_dir} for future use? [Y/n]: ").lower() + user_consent_to_save = response in ('y', 'yes', '') + + if user_consent_to_save: + try: + cache_dir.mkdir(parents=True, exist_ok=True) + colour.write_LUT(colour.LUT3D(table=generated_exposure_lut_table, name=f"{base_filename}_exposure"), str(exposure_lut_path)) + print(f"Saved exposure LUT to {exposure_lut_path}") + except Exception as e: + print(f"Error saving exposure LUT to {exposure_lut_path}: {e}", file=sys.stderr) + + # Ensure exposure_lut is now populated. If not, something went wrong. + if exposure_lut is None: + # This case should ideally not be reached if create_exposure_lut always returns a valid LUT or raises an error. + print("Critical error: Exposure LUT could not be loaded or generated. Exiting.", file=sys.stderr) + sys.exit(1) + + # For density LUTs, we need to know the typical range of log exposure values + # We can estimate this from the H&D curve data itself. + log_e_points = [p.d for p in datasheet.properties.curves.hd if p.d is not None] + log_exposure_range = (log_e_points[0], log_e_points[-1]) + + density_lut_naive: Optional[np.ndarray] = None # Initialize naive density LUT + if not args.no_cache and naive_density_lut_path.exists(): + try: + print(f"Loading naive density LUT from cache: {naive_density_lut_path}") + loaded_obj = colour.read_LUT(str(naive_density_lut_path)) + if isinstance(loaded_obj, colour.LUT3D): + density_lut_naive = loaded_obj.table + print(f"Successfully loaded naive density LUT from {naive_density_lut_path}") + else: + # Log a warning and fall through to regeneration + print(f"Warning: Cached naive density LUT {naive_density_lut_path} is not a LUT3D (type: {type(loaded_obj)}). Will regenerate.") + except Exception as e: + # Log a warning and fall through to regeneration + print(f"Warning: Error loading naive density LUT from {naive_density_lut_path}: {e}. Will regenerate.") + + if density_lut_naive is None: + print("Generating new naive density LUT...") + density_lut_naive = create_density_lut( + datasheet.processing, + datasheet.properties.curves.hd, + middle_gray_logE, + log_exposure_range, + size=LUT_SIZE, + ) + if not args.no_cache: + if user_consent_to_save is None: + response = input(f"Save generated naive density LUT to {cache_dir} for future use? [Y/n]: ").lower() + user_consent_to_save = response in ('y', 'yes', '') + if user_consent_to_save: + try: + cache_dir.mkdir(parents=True, exist_ok=True) + colour.write_LUT( + colour.LUT3D(table=density_lut_naive, name=f"{base_filename}_naive_density"), + str(naive_density_lut_path), + ) + print(f"Saved naive density LUT to {naive_density_lut_path}") + except Exception as e: + print(f"Error saving naive density LUT to {naive_density_lut_path}: {e}", file=sys.stderr) + + if density_lut_naive is None: + # This case should ideally not be reached if create_density_lut always returns a valid LUT or raises an error. + print("Critical error: Naive density LUT could not be loaded or generated. Exiting.", file=sys.stderr) + sys.exit(1) + + + inhibitor_matrix = compute_inhibitor_matrix( + datasheet.properties.couplers.dir_amount_rgb, + datasheet.properties.couplers.dir_diffusion_interlayer, + ) + + if not args.no_cache and uncoupled_density_lut_path.exists(): + try: + print(f"Loading uncoupled density LUT from cache: {uncoupled_density_lut_path}") + loaded_obj = colour.read_LUT(str(uncoupled_density_lut_path)) + if isinstance(loaded_obj, colour.LUT3D): + density_lut_uncoupled = loaded_obj.table + print(f"Successfully loaded uncoupled density LUT from {uncoupled_density_lut_path}") + else: + # Log a warning and fall through to regeneration + print(f"Warning: Cached uncoupled density LUT {uncoupled_density_lut_path} is not a LUT3D (type: {type(loaded_obj)}). Will regenerate.") + except Exception as e: + # Log a warning and fall through to regeneration + print(f"Warning: Error loading uncoupled density LUT from {uncoupled_density_lut_path}: {e}. Will regenerate.") + else: + uncoupled_hd_curve = compute_uncoupled_hd_curves( + datasheet.properties.curves.hd, inhibitor_matrix + ) + density_lut_uncoupled = create_density_lut( + datasheet.processing, + uncoupled_hd_curve, + middle_gray_logE, + log_exposure_range, + size=LUT_SIZE, + ) + if not args.no_cache: + if user_consent_to_save is None: + response = input(f"Save generated uncoupled density LUT to {cache_dir} for future use? [Y/n]: ").lower() + user_consent_to_save = response in ('y', 'yes', '') + if user_consent_to_save: + try: + cache_dir.mkdir(parents=True, exist_ok=True) + colour.write_LUT( + colour.LUT3D(table=density_lut_uncoupled, name=f"{base_filename}_uncoupled_density"), + str(uncoupled_density_lut_path), + ) + print(f"Saved uncoupled density LUT to {uncoupled_density_lut_path}") + except Exception as e: + print(f"Error saving uncoupled density LUT to {uncoupled_density_lut_path}: {e}", file=sys.stderr) + + if density_lut_uncoupled is None: + # This case should ideally not be reached if create_density_lut always returns a valid LUT or raises an error. + print("Critical error: Uncoupled density LUT could not be loaded or generated. Exiting.", file=sys.stderr) + sys.exit(1) + + # --- Load and Prepare 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") + if any( + args.input_image.lower().endswith(ext) for ext in [".dng", ".raw", ".arw"] ): - print("Detected Camera RAW file, reading raw image data...") with rawpy.imread(args.input_image) as raw: + if args.force_d65: + wb_params = {'user_wb': (1.0, 1.0, 1.0, 1.0)} # D65 white balance 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 + demosaic_algorithm=rawpy.DemosaicAlgorithm.AHD, # type: ignore + output_bps=16, + gamma=(1.0, 1.0), + no_auto_bright=True, + output_color=rawpy.ColorSpace.sRGB, # type: ignore + **(wb_params if args.force_d65 else {'user_camera_wb': True}), ) - # 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] + image_linear = image_raw.astype(np.float64) / 65535.0 + print("RAW file processed to linear floating point directly.") 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) + image_float = image_raw.astype(np.float64) / ( + 65535.0 if image_raw.dtype == np.uint16 else 255.0 + ) + image_linear = ( + colour.cctf_decoding(image_float, function="sRGB") + if image_raw.dtype in (np.uint8, np.uint16) + else image_float + ) + if args.force_d65: + image_linear = white_balance_to_d65(image_linear) 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.") + if image_linear.shape[-1] > 3: 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] + image_linear, image_width_px = np.maximum(image_linear, 0.0), image_linear.shape[1] # --- Pipeline Steps --- + ### PERFORMANCE ### + # 1. Convert Linear RGB to Log Exposure using the LUT + print("Applying exposure LUT...") + log_exposure_rgb = apply_3d_lut(image_linear, exposure_lut) - # 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) - - + # 2. Apply DIR Coupler Simulation 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 + # 2a. Calculate "naive" density using its LUT + naive_density_rgb = apply_3d_lut( + (log_exposure_rgb - log_exposure_range[0]) + / (log_exposure_range[1] - log_exposure_range[0]), + density_lut_naive, + ) + # naive_density_rgb = apply_hd_curves(log_exposure_rgb, datasheet.processing, datasheet.properties.curves.hd, middle_gray_logE) # Apply H&D curves to naive density + + # 2b. Use this density to modify 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 + 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. Apply the *uncoupled* density LUT to the *modified* exposure + print("Applying uncoupled density LUT...") + # Normalize the modified log exposure to the [0,1] range for LUT lookup + norm_log_exposure = (modified_log_exposure_rgb - log_exposure_range[0]) / ( + log_exposure_range[1] - log_exposure_range[0] ) + density_rgb = apply_3d_lut(np.clip(norm_log_exposure, 0, 1), density_lut_uncoupled) - # 3. Convert Density back to Linear Transmittance - # Higher density means lower transmittance + # (Rest of the pipeline is unchanged) 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) + linear_transmittance = np.clip(10.0 ** (-density_rgb), 0.0, 1.0) 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.couplers, 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) + print("Applying saturation adjustment...") + linear_post_saturation = apply_saturation_rgb( + linear_post_spatial, datasheet.properties.couplers.saturation_amount + ) + print("Converting to output format and saving...") + output_image = ( + np.clip(linear_post_saturation, 0.0, 1.0) + * (65535.0 if args.output_image.lower().endswith((".tiff", ".tif")) else 255.0) + ).astype( + np.uint16 if args.output_image.lower().endswith((".tiff", ".tif")) else np.uint8 + ) + iio.imwrite(args.output_image, output_image) + print("Done.") if __name__ == "__main__": diff --git a/pyproject.toml b/pyproject.toml index d10842e..9d5c67d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,7 @@ dependencies = [ "jupyterlab>=4.4.3", "numpy>=2.2.6", "pillow>=11.2.1", + "pyfftw>=0.15.0", "rawpy>=0.25.0", "scipy>=1.15.3", "warp-lang>=1.7.2", diff --git a/sim_data/ektar_100.json b/sim_data/ektar_100.json index a2cd73d..32427b8 100644 --- a/sim_data/ektar_100.json +++ b/sim_data/ektar_100.json @@ -20,7 +20,7 @@ "properties": { "calibration": { "iso": 100, - "middle_gray_logh": -0.84 + "middle_gray_logE": -0.84 }, "halation": { "strength": { diff --git a/sim_data/portra_400.json b/sim_data/portra_400.json index 9e50005..d34cc13 100644 --- a/sim_data/portra_400.json +++ b/sim_data/portra_400.json @@ -20,7 +20,7 @@ "properties": { "calibration": { "iso": 400, - "middle_gray_logh": -1.44 + "middle_gray_logE": -1.44 }, "halation": { "strength": { diff --git a/test_images/RAW/07.DNG.cfg b/test_images/RAW/07.DNG.cfg new file mode 100644 index 0000000..e7b6450 --- /dev/null +++ b/test_images/RAW/07.DNG.cfg @@ -0,0 +1,84 @@ +frames:1 +fps:0 +module:i-raw:main:43:400 +module:denoise:01:218:400 +module:hilite:01:393:400 +module:demosaic:01:568:400 +module:colour:01:931:400 +module:filmcurv:01:1094:400 +module:llap:01:1269:400 +module:hist:01:1418:634 +module:zones:01:219:800 +module:crop:01:743:400 +module:lens:01:568:800 +module:pick:01:43:800 +module:display:hist:1600:634 +module:display:main:1545:400 +connect:i-raw:main:output:denoise:01:input +connect:denoise:01:output:hilite:01:input +connect:hilite:01:output:demosaic:01:input +connect:crop:01:output:colour:01:input +connect:colour:01:output:filmcurv:01:input +connect:filmcurv:01:output:llap:01:input +connect:llap:01:output:hist:01:input +connect:demosaic:01:output:crop:01:input +connect:hist:01:output:display:hist:input +connect:llap:01:output:display:main:input +param:i-raw:main:filename:07.DNG +param:i-raw:main:noise a:0 +param:i-raw:main:noise b:0 +param:i-raw:main:startid:0 +param:denoise:01:strength:0 +param:denoise:01:luma:0.6 +param:denoise:01:detail:1 +param:denoise:01:gainmap:1 +param:hilite:01:white:0.985 +param:hilite:01:desat:0.3 +param:hilite:01:soft:0.6 +param:demosaic:01:colour:0 +param:demosaic:01:method:0 +param:colour:01:exposure:0 +param:colour:01:sat:1 +param:colour:01:picked:0 +param:colour:01:matrix:1 +param:colour:01:gamut:0 +param:colour:01:clip:0 +param:colour:01:clipmax:1 +param:colour:01:temp:6504 +param:colour:01:white:1.19936:1:1.04374:0 +param:colour:01:mat:-1:0:0:0:0:0:0:0:0 +param:colour:01:mode:0 +param:colour:01:cnt:4 +param:colour:01:rbmap:0.3333:0.3333:0.3333:0.3333:0.3333:0.3333:0.5:0.25:0.25:0.5:0.25:0.25:0.25:0.5:0.25:0.25:0.5:0.25:0.25:0.25:0.5:0.25:0.25:0.5:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 +param:colour:01:import:01 +param:filmcurv:01:light:3 +param:filmcurv:01:contrast:1.2 +param:filmcurv:01:bias:0 +param:filmcurv:01:colour:3 +param:llap:01:sigma:0.12 +param:llap:01:shadows:1 +param:llap:01:hilights:1 +param:llap:01:clarity:0.2 +param:zones:01:radius:0.01 +param:zones:01:epsilon:0.06 +param:zones:01:gamma:1 +param:zones:01:nzones:3 +param:zones:01:zone:0:0:0:0:0:0:0 +param:crop:01:perspect:0.25:0.25:0.75:0.25:0.75:0.75:0.25:0.75 +param:crop:01:crop:0.000315126:0.999685:0.000473485:0.999527 +param:crop:01:rotate:0 +param:lens:01:center:0:0 +param:lens:01:scale:1:1 +param:lens:01:squish0:0 +param:lens:01:squish1:0 +param:lens:01:ca red:1 +param:lens:01:ca blue:1 +param:pick:01:nspots:0 +param:pick:01:pad:0:0:0 +param:pick:01:spots:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 +param:pick:01:picked:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 +param:pick:01:ref:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 +param:pick:01:show:0 +param:pick:01:grab:0 +param:pick:01:de76:0:0:0 +param:pick:01:freeze:0 diff --git a/test_images/v1.1output/filmcolor/09.tiff.jpg b/test_images/v1.1output/filmcolor/09.tiff.jpg new file mode 100644 index 0000000..b4f2519 --- /dev/null +++ b/test_images/v1.1output/filmcolor/09.tiff.jpg @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c3df25694eb483f63a2d865e7a1b4838e0bb21a85a23bce4ac99d6bc1d5aa13d +size 7229760 diff --git a/test_images/v1.1output/filmscan/09.tiff.jpg b/test_images/v1.1output/filmscan/09.tiff.jpg new file mode 100644 index 0000000..edb3eb1 --- /dev/null +++ b/test_images/v1.1output/filmscan/09.tiff.jpg @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6fca48d71a2b9c4b161cadfd1e24d070863168de04b0eba9597a8e983e489592 +size 9542730 diff --git a/test_images/v1.3output/filmcolor/diff.png b/test_images/v1.3output/filmcolor/diff.png new file mode 100644 index 0000000..23fd68a Binary files /dev/null and b/test_images/v1.3output/filmcolor/diff.png differ diff --git a/test_images/v1.3output/filmcolor/diff2.png b/test_images/v1.3output/filmcolor/diff2.png new file mode 100644 index 0000000..9664a6a Binary files /dev/null and b/test_images/v1.3output/filmcolor/diff2.png differ diff --git a/test_images/v1.3output/filmcolor/diff3.png b/test_images/v1.3output/filmcolor/diff3.png new file mode 100644 index 0000000..a8cf651 Binary files /dev/null and b/test_images/v1.3output/filmcolor/diff3.png differ diff --git a/test_images/v1.3output/filmcolor/diff4.png b/test_images/v1.3output/filmcolor/diff4.png new file mode 100644 index 0000000..ec00514 Binary files /dev/null and b/test_images/v1.3output/filmcolor/diff4.png differ diff --git a/test_images/v1.3output/filmcolor/diff5.png b/test_images/v1.3output/filmcolor/diff5.png new file mode 100644 index 0000000..b5b12b2 Binary files /dev/null and b/test_images/v1.3output/filmcolor/diff5.png differ diff --git a/test_images/v1.3output/filmcolor/ektar_diff_wb.png b/test_images/v1.3output/filmcolor/ektar_diff_wb.png new file mode 100644 index 0000000..8e60da7 Binary files /dev/null and b/test_images/v1.3output/filmcolor/ektar_diff_wb.png differ diff --git a/uv.lock b/uv.lock index fdde432..ba1f7ce 100644 --- a/uv.lock +++ b/uv.lock @@ -285,6 +285,7 @@ dependencies = [ { name = "jupyterlab" }, { name = "numpy" }, { name = "pillow" }, + { name = "pyfftw" }, { name = "rawpy" }, { name = "scipy" }, { name = "warp-lang" }, @@ -298,6 +299,7 @@ requires-dist = [ { name = "jupyterlab", specifier = ">=4.4.3" }, { name = "numpy", specifier = ">=2.2.6" }, { name = "pillow", specifier = ">=11.2.1" }, + { name = "pyfftw", specifier = ">=0.15.0" }, { name = "rawpy", specifier = ">=0.25.0" }, { name = "scipy", specifier = ">=1.15.3" }, { name = "warp-lang", specifier = ">=1.7.2" }, @@ -1056,6 +1058,24 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/13/a3/a812df4e2dd5696d1f351d58b8fe16a405b234ad2886a0dab9183fb78109/pycparser-2.22-py3-none-any.whl", hash = "sha256:c3702b6d3dd8c7abc1afa565d7e63d53a1d0bd86cdc24edd75470f4de499cfcc", size = 117552 }, ] +[[package]] +name = "pyfftw" +version = "0.15.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy" }, + { name = "setuptools" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/4b/3f/ee1bc44b080fc1e81d293cd07bed563d254bc1997d63a3b8053804a87dfd/pyfftw-0.15.0.tar.gz", hash = "sha256:2f16b9854a40c8fdd10aa5803b24ddc6ab49f9cd559dbd7f07e7d61aa205c1ca", size = 164003 } +wheels = [ + { url = "https://files.pythonhosted.org/packages/a1/0d/90821a720638cbef79319463e3ce5d98f24153cee53316d4ecc1d8c0400d/pyFFTW-0.15.0-cp313-cp313-macosx_12_0_x86_64.whl", hash = "sha256:cfbeb106877db1b6bf527735647b861c86ac846ff47671d0d855a9be2de83368", size = 2832851 }, + { url = "https://files.pythonhosted.org/packages/12/b7/e625f0cdb2bd65aa43ee35ce8cac1b0abc6c12871cd87fa83e404d82fd1d/pyFFTW-0.15.0-cp313-cp313-macosx_14_0_arm64.whl", hash = "sha256:6d8913036a48ebcc9e3e1a315a6607e5cc31af4aee395aae180ff644c4658bbb", size = 1269651 }, + { url = "https://files.pythonhosted.org/packages/0f/66/d12a9629908008ed97446c44ca2177c8f59e570a445da4136fdbb9e8db1c/pyFFTW-0.15.0-cp313-cp313-manylinux_2_28_aarch64.whl", hash = "sha256:05324c8c092aa43868554daf6ddd7fd23eee18cd284ab68e749c62427c662737", size = 2496935 }, + { url = "https://files.pythonhosted.org/packages/ad/c4/adcf010b151564401e3d401602b0e0cacae1e7cb816f0207397c3ddec626/pyFFTW-0.15.0-cp313-cp313-manylinux_2_28_x86_64.whl", hash = "sha256:474f4a0ccfddfdfbe3e4a610b1c8f5968b58edf44434bf9d08ea52108a72df54", size = 3096170 }, + { url = "https://files.pythonhosted.org/packages/b2/96/b1eaad07d5eaf026e2ed0c6f2afdf825628706162071d2665b3824727346/pyFFTW-0.15.0-cp313-cp313-win32.whl", hash = "sha256:40f8f3341546264178a8d9e8736e91554884595a683b71f8db8399907330f47d", size = 2229713 }, + { url = "https://files.pythonhosted.org/packages/e9/23/d06a3e5b549f7537ada2b1fc2e5f7c05c38df70c843d463bc469812b1ebd/pyFFTW-0.15.0-cp313-cp313-win_amd64.whl", hash = "sha256:bf04458c5d2fbfe5da270a4b667c59c5792cbe39300379d886acd8ea97fc55cc", size = 2640821 }, +] + [[package]] name = "pygments" version = "2.19.1"