From 0b0b32356c1c83b7036c204b372fd78d54eecc64 Mon Sep 17 00:00:00 2001 From: Tanishq Dubey Date: Sat, 7 Jun 2025 19:29:33 -0400 Subject: [PATCH] better dir sim --- filmcolor | 159 ++++++++++++++++++++++++++++++--------- sim_data/ektar_100.json | 6 +- sim_data/portra_400.json | 6 +- 3 files changed, 131 insertions(+), 40 deletions(-) diff --git a/filmcolor b/filmcolor index 87bfbd2..2c5df1b 100755 --- a/filmcolor +++ b/filmcolor @@ -6,7 +6,7 @@ # "Pillow", # "imageio", # "rawpy", -# "colour", +# "colour-science", # ] # /// @@ -26,7 +26,7 @@ import csv import numpy as np import imageio.v3 as iio from scipy.interpolate import interp1d -from scipy.ndimage import gaussian_filter +from scipy.ndimage import gaussian_filter, gaussian_filter1d import rawpy import os import sys @@ -110,15 +110,16 @@ class Processing: class Couplers: - amount: float - diffusion_um: float + saturation_amount: float + dir_amount_rgb: List[float] + dir_diffusion_um: float + dir_diffusion_interlayer: float - def __init__(self, amount: float, diffusion_um: float) -> None: - self.amount = amount - self.diffusion_um = diffusion_um - - def __repr__(self) -> str: - return f"Couplers(amount={self.amount:.3f}, diffusion_um={self.diffusion_um:.1f})" + 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: @@ -309,8 +310,10 @@ def parse_datasheet_json(json_filepath) -> FilmDatasheet | None: ), ) couplers = Couplers( - amount=data["properties"]["couplers"]["amount"], - diffusion_um=data["properties"]["couplers"]["diffusion_um"], + 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"] @@ -359,6 +362,87 @@ def um_to_pixels(sigma_um, image_width_px, film_format_mm): 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, @@ -465,32 +549,19 @@ def apply_saturation_rgb(image_linear, saturation_factor): 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) - total_diffusion_um = np.sqrt( - couplerData.diffusion_um**2 + interlayerData.diffusion_um**2 - ) - - if total_diffusion_um > EPSILON: - sigma_pixels_diffusion = um_to_pixels( - total_diffusion_um, image_width_px, film_format_mm - ) + 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 diffusion blur: sigma={sigma_pixels_diffusion:.2f} pixels ({total_diffusion_um:.1f} um)" - ) - # Apply blur to the linear image data - image = gaussian_filter( - image, - sigma=[sigma_pixels_diffusion, sigma_pixels_diffusion, 0], - mode="nearest", - ) # Blur R, G, B independently - image = np.clip(image, 0.0, 1.0) # Keep values in range + 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 @@ -643,6 +714,12 @@ 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.") + import pprint pprint.pp(datasheet) @@ -740,12 +817,25 @@ def main(): 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( - log_exposure_rgb, + modified_log_exposure_rgb, datasheet.processing, datasheet.properties.curves.hd, middle_gray_logE, @@ -767,18 +857,15 @@ def main(): 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) - print("Applying saturation adjustment...") - coupler_amount = datasheet.properties.couplers.amount # Assuming coupler_amount directly scales saturation factor. # Values > 1 increase saturation, < 1 decrease. - linear_post_saturation = apply_saturation_rgb(linear_post_spatial, coupler_amount) + linear_post_saturation = apply_saturation_rgb(linear_post_spatial, datasheet.properties.couplers.saturation_amount) # --- Final Output Conversion --- print("Converting to output format...") diff --git a/sim_data/ektar_100.json b/sim_data/ektar_100.json index 3e5d227..a2cd73d 100644 --- a/sim_data/ektar_100.json +++ b/sim_data/ektar_100.json @@ -35,8 +35,10 @@ } }, "couplers": { - "amount": 1.0, - "diffusion_um": 5.0 + "saturation_amount": 1.1, + "dir_amount_rgb": [0.7, 0.9, 0.5], + "dir_diffusion_um": 11.0, + "dir_diffusion_interlayer": 1.9 }, "interlayer": { "diffusion_um": 2.1 diff --git a/sim_data/portra_400.json b/sim_data/portra_400.json index d99a67f..9e50005 100644 --- a/sim_data/portra_400.json +++ b/sim_data/portra_400.json @@ -35,8 +35,10 @@ } }, "couplers": { - "amount": 1.0, - "diffusion_um": 5.0 + "saturation_amount": 1.0, + "dir_amount_rgb": [0.7, 0.9, 0.5], + "dir_diffusion_um": 15.0, + "dir_diffusion_interlayer": 1.5 }, "interlayer": { "diffusion_um": 2.1