Raw-rs: add post-processing steps (#1923)

* add convert_to_rgb step

* add code to generate gamma correction curve

* add gamma correction step

* fix clippy warnings and cargo fmt

* remove unnecessary dependencies

* Code review 1

* Code review 2

* fix the order of operations

* Code review 3

---------

Co-authored-by: Keavon Chambers <keavon@keavon.com>

Cargo.lock

@@ -743,12 +743,6 @@ version = "3.16.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "79296716171880943b8470b5f8d03aa55eb2e645a4874bdbb28adb49162e012c"
 
-[[package]]
-name = "by_address"
-version = "1.2.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "64fa3c856b712db6612c019f14756e64e4bcea13337a6b33b696333a9eaa2d06"
-
 [[package]]
 name = "bytemuck"
 version = "1.16.1"
@@ -1695,12 +1689,6 @@ dependencies = [
  "pin-project-lite",
 ]
 
-[[package]]
-name = "fast-srgb8"
-version = "1.0.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "dd2e7510819d6fbf51a5545c8f922716ecfb14df168a3242f7d33e0239efe6a1"
-
 [[package]]
 name = "fastnoise-lite"
 version = "1.1.1"
@@ -4351,30 +4339,6 @@ dependencies = [
  "ttf-parser 0.24.0",
 ]
 
-[[package]]
-name = "palette"
-version = "0.7.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4cbf71184cc5ecc2e4e1baccdb21026c20e5fc3dcf63028a086131b3ab00b6e6"
-dependencies = [
- "approx",
- "fast-srgb8",
- "palette_derive",
- "phf 0.11.2",
-]
-
-[[package]]
-name = "palette_derive"
-version = "0.7.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f5030daf005bface118c096f510ffb781fc28f9ab6a32ab224d8631be6851d30"
-dependencies = [
- "by_address",
- "proc-macro2",
- "quote",
- "syn 2.0.72",
-]
-
 [[package]]
 name = "pango"
 version = "0.15.10"
@@ -5003,7 +4967,6 @@ dependencies = [
  "image 0.25.2",
  "libraw-rs",
  "num_enum 0.7.2",
- "palette",
  "reqwest 0.12.5",
  "tag-derive",
  "thiserror",

libraries/raw-rs/Cargo.toml

@@ -13,7 +13,7 @@ repository = "https://github.com/GraphiteEditor/Graphite/tree/master/libraries/r
 documentation = "https://docs.rs/raw-rs"
 
 [features]
-raw-rs-tests = ["dep:image", "dep:libraw-rs", "dep:palette", "dep:reqwest"]
+raw-rs-tests = ["dep:image", "dep:libraw-rs", "dep:reqwest"]
 
 [dependencies]
 # Local dependencies
@@ -33,4 +33,3 @@ reqwest = { workspace = true, optional = true }
 
 # Optional dependencies (should be dev dependencies, but Cargo currently doesn't allow optional dev dependencies)
 libraw-rs = { version = "0.0.4", optional = true }
-palette = { version = "0.7.6", optional = true }

libraries/raw-rs/src/decoder/arw1.rs

@@ -26,7 +26,10 @@ pub fn decode_a100<R: Read + Seek>(ifd: Ifd, file: &mut TiffRead<R>) -> RawImage
 		#[allow(unreachable_code)]
 		maximum: (1 << 12) - 1,
 		black: SubtractBlack::None,
-		camera_to_xyz: None,
+		camera_model: None,
+		white_balance_multiplier: None,
+		camera_to_rgb: None,
+		rgb_to_camera: None,
 	}
 }
 

libraries/raw-rs/src/decoder/arw2.rs

@@ -49,7 +49,10 @@ pub fn decode<R: Read + Seek>(ifd: Ifd, file: &mut TiffRead<R>) -> RawImage {
 		cfa_pattern: ifd.cfa_pattern.try_into().unwrap(),
 		maximum: (1 << 14) - 1,
 		black: SubtractBlack::CfaGrid([512, 512, 512, 512]), // TODO: Find the correct way to do this
-		camera_to_xyz: None,
+		camera_model: None,
+		white_balance_multiplier: None,
+		camera_to_rgb: None,
+		rgb_to_camera: None,
 	}
 }
 

libraries/raw-rs/src/decoder/uncompressed.rs

@@ -57,6 +57,9 @@ pub fn decode<R: Read + Seek>(ifd: Ifd, file: &mut TiffRead<R>) -> RawImage {
 		cfa_pattern: ifd.cfa_pattern.try_into().unwrap(),
 		maximum: if bits_per_sample == 16 { u16::MAX } else { (1 << bits_per_sample) - 1 },
 		black: SubtractBlack::CfaGrid(ifd.black_level),
-		camera_to_xyz: None,
+		camera_model: None,
+		white_balance_multiplier: None,
+		camera_to_rgb: None,
+		rgb_to_camera: None,
 	}
 }

libraries/raw-rs/src/demosaicing/linear_demosaicing.rs

@@ -73,5 +73,7 @@ fn linear_demosaic_rggb(mut raw_image: RawImage) -> Image<u16> {
 		data: raw_image.data,
 		width: raw_image.width,
 		height: raw_image.height,
+		rgb_to_camera: raw_image.rgb_to_camera,
+		histogram: None,
 	}
 }

libraries/raw-rs/src/lib.rs

@@ -1,10 +1,11 @@
 pub mod decoder;
 pub mod demosaicing;
 pub mod metadata;
+pub mod postprocessing;
 pub mod preprocessing;
 pub mod tiff;
 
-use crate::preprocessing::camera_data::camera_to_xyz;
+use crate::metadata::identify::CameraModel;
 
 use tag_derive::Tag;
 use tiff::file::TiffRead;
@@ -27,14 +28,21 @@ pub struct RawImage {
 	pub cfa_pattern: [u8; 4],
 	pub maximum: u16,
 	pub black: SubtractBlack,
-	pub camera_to_xyz: Option<[f64; 9]>,
+	pub camera_model: Option<CameraModel>,
+	pub white_balance_multiplier: Option<[f64; 3]>,
+	pub camera_to_rgb: Option<[[f64; 3]; 3]>,
+	pub rgb_to_camera: Option<[[f64; 3]; 3]>,
 }
 
 pub struct Image<T> {
 	pub data: Vec<T>,
 	pub width: usize,
 	pub height: usize,
+	/// We can assume this will be 3 for all non-obscure, modern cameras.
+	/// See <https://github.com/GraphiteEditor/Graphite/pull/1923#discussion_r1725070342> for more information.
 	pub channels: u8,
+	pub rgb_to_camera: Option<[[f64; 3]; 3]>,
+	pub histogram: Option<[[usize; 0x2000]; 3]>,
 }
 
 #[allow(dead_code)]
@@ -68,27 +76,32 @@ pub fn decode<R: Read + Seek>(reader: &mut R) -> Result<RawImage, DecoderError>
 		}
 	};
 
-	raw_image.camera_to_xyz = camera_to_xyz(&camera_model);
+	raw_image.camera_model = Some(camera_model);
 
 	Ok(raw_image)
 }
 
 pub fn process_8bit(raw_image: RawImage) -> Image<u8> {
-	let raw_image = crate::preprocessing::subtract_black::subtract_black(raw_image);
-	let raw_image = crate::preprocessing::raw_to_image::raw_to_image(raw_image);
-	let raw_image = crate::preprocessing::scale_colors::scale_colors(raw_image);
-	let image = crate::demosaicing::linear_demosaicing::linear_demosaic(raw_image);
+	let image = process_16bit(raw_image);
 
 	Image {
 		channels: image.channels,
 		data: image.data.iter().map(|x| (x >> 8) as u8).collect(),
 		width: image.width,
 		height: image.height,
+		rgb_to_camera: image.rgb_to_camera,
+		histogram: image.histogram,
 	}
 }
 
-pub fn process_16bit(_image: RawImage) -> Image<u16> {
-	todo!()
+pub fn process_16bit(raw_image: RawImage) -> Image<u16> {
+	let raw_image = crate::preprocessing::camera_data::calculate_conversion_matrices(raw_image);
+	let raw_image = crate::preprocessing::subtract_black::subtract_black(raw_image);
+	let raw_image = crate::preprocessing::raw_to_image::raw_to_image(raw_image);
+	let raw_image = crate::preprocessing::scale_colors::scale_colors(raw_image);
+	let image = crate::demosaicing::linear_demosaicing::linear_demosaic(raw_image);
+	let image = crate::postprocessing::convert_to_rgb::convert_to_rgb(image);
+	crate::postprocessing::gamma_correction::gamma_correction(image)
 }
 
 #[derive(Error, Debug)]

libraries/raw-rs/src/postprocessing/convert_to_rgb.rs

@@ -0,0 +1,40 @@
+use crate::Image;
+
+const CHANNELS_IN_RGB: usize = 3;
+
+pub fn convert_to_rgb(mut image: Image<u16>) -> Image<u16> {
+	let Some(rgb_to_camera) = image.rgb_to_camera else { return image };
+
+	// Rarely this might be 4 instead of 3 if an obscure Bayer filter is used, such as RGBE or CYGM, instead of the typical RGGB.
+	// See: <https://github.com/GraphiteEditor/Graphite/pull/1923#discussion_r1725070342>.
+	let channels = image.channels as usize;
+	let mut data = Vec::with_capacity(CHANNELS_IN_RGB * image.width * image.height);
+	let mut histogram = [[0; 0x2000]; CHANNELS_IN_RGB];
+
+	for i in 0..(image.height * image.width) {
+		let start = i * channels;
+		let end = start + channels;
+		let input_pixel = &mut image.data[start..end];
+
+		let mut output_pixel = [0.; CHANNELS_IN_RGB];
+		for (channel, &value) in input_pixel.iter().enumerate() {
+			output_pixel[0] += rgb_to_camera[0][channel] * value as f64;
+			output_pixel[1] += rgb_to_camera[1][channel] * value as f64;
+			output_pixel[2] += rgb_to_camera[2][channel] * value as f64;
+		}
+
+		for (output_pixel_channel, histogram_channel) in output_pixel.iter().zip(histogram.iter_mut()) {
+			let final_sum = (*output_pixel_channel as u16).clamp(0, u16::MAX);
+
+			histogram_channel[final_sum as usize >> CHANNELS_IN_RGB] += 1;
+
+			data.push(final_sum);
+		}
+	}
+
+	image.data = data;
+	image.histogram = Some(histogram);
+	image.channels = CHANNELS_IN_RGB as u8;
+
+	image
+}

libraries/raw-rs/src/postprocessing/gamma_correction.rs

@@ -0,0 +1,89 @@
+use crate::Image;
+use std::f64::consts::E;
+
+pub fn gamma_correction(mut image: Image<u16>) -> Image<u16> {
+	let Some(histogram) = image.histogram else { return image };
+
+	let percentage = image.width * image.height;
+
+	let mut white = 0;
+	for channel_histogram in histogram {
+		let mut total = 0;
+		for i in (0x20..0x2000).rev() {
+			total += channel_histogram[i] as u64;
+
+			if total * 100 > percentage as u64 {
+				white = white.max(i);
+				break;
+			}
+		}
+	}
+
+	let curve = generate_gamma_curve(0.45, 4.5, (white << 3) as f64);
+
+	for value in image.data.iter_mut() {
+		*value = curve[*value as usize];
+	}
+	image.histogram = None;
+
+	image
+}
+
+/// `max_intensity` must be non-zero.
+fn generate_gamma_curve(power: f64, threshold: f64, max_intensity: f64) -> Vec<u16> {
+	debug_assert!(max_intensity != 0.);
+
+	let (mut bound_start, mut bound_end) = if threshold >= 1. { (0., 1.) } else { (1., 0.) };
+
+	let mut transition_point = 0.;
+	let mut transition_ratio = 0.;
+	let mut curve_adjustment = 0.;
+
+	if threshold != 0. && (threshold - 1.) * (power - 1.) <= 0. {
+		for _ in 0..48 {
+			transition_point = (bound_start + bound_end) / 2.;
+
+			if power != 0. {
+				let temp_transition_ratio = transition_point / threshold;
+				let exponential_power = temp_transition_ratio.powf(-power);
+				let normalized_exponential_power = (exponential_power - 1.) / power;
+				let comparison_result = normalized_exponential_power - (1. / transition_point);
+
+				let bound_to_update = if comparison_result > -1. { &mut bound_end } else { &mut bound_start };
+				*bound_to_update = transition_point;
+			} else {
+				let adjusted_transition_point = E.powf(1. - 1. / transition_point);
+				let transition_point_ratio = transition_point / adjusted_transition_point;
+
+				let bound_to_update = if transition_point_ratio < threshold { &mut bound_end } else { &mut bound_start };
+				*bound_to_update = transition_point;
+			}
+		}
+
+		transition_ratio = transition_point / threshold;
+
+		if power != 0. {
+			curve_adjustment = transition_point * ((1. / power) - 1.);
+		}
+	}
+
+	let mut curve = vec![0xffff; 0x1_0000];
+	let length = curve.len() as f64;
+
+	for (i, entry) in curve.iter_mut().enumerate() {
+		let ratio = (i as f64) / max_intensity;
+		if ratio < 1. {
+			let altered_ratio = if ratio < transition_ratio {
+				ratio * threshold
+			} else if power != 0. {
+				ratio.powf(power) * (1. + curve_adjustment) - curve_adjustment
+			} else {
+				ratio.ln() * transition_point + 1.
+			};
+
+			*entry = (length * altered_ratio) as u16;
+		}
+	}
+
+	curve
+}

libraries/raw-rs/src/postprocessing/mod.rs

@@ -0,0 +1,2 @@
+pub mod convert_to_rgb;
+pub mod gamma_correction;

libraries/raw-rs/src/preprocessing/camera_data.rs

@@ -1,4 +1,4 @@
-use crate::metadata::identify::CameraModel;
+use crate::RawImage;
 use build_camera_data::build_camera_data;
 
 pub struct CameraData {
@@ -17,10 +17,97 @@ impl CameraData {
 
 const CAMERA_DATA: [(&str, CameraData); 40] = build_camera_data!();
 
-pub fn camera_to_xyz(camera_model: &CameraModel) -> Option<[f64; 9]> {
+const XYZ_TO_RGB: [[f64; 3]; 3] = [
+	// Matrix:
+	[0.412453, 0.357580, 0.180423],
+	[0.212671, 0.715160, 0.072169],
+	[0.019334, 0.119193, 0.950227],
+];
+
+pub fn calculate_conversion_matrices(mut raw_image: RawImage) -> RawImage {
+	let Some(ref camera_model) = raw_image.camera_model else { return raw_image };
 	let camera_name_needle = camera_model.make.to_owned() + " " + &camera_model.model;
-	CAMERA_DATA
+
+	let camera_to_xyz = CAMERA_DATA
 		.iter()
 		.find(|(camera_name_haystack, _)| camera_name_needle == *camera_name_haystack)
-		.map(|(_, data)| data.camera_to_xyz.map(|x| (x as f64) / 10_000.))
+		.map(|(_, data)| data.camera_to_xyz.map(|x| (x as f64) / 10_000.));
+	let Some(camera_to_xyz) = camera_to_xyz else { return raw_image };
+
+	let mut camera_to_rgb = [[0.; 3]; 3];
+	for i in 0..3 {
+		for j in 0..3 {
+			for k in 0..3 {
+				camera_to_rgb[i][j] += camera_to_xyz[i * 3 + k] * XYZ_TO_RGB[k][j];
+			}
+		}
+	}
+
+	let white_balance_multiplier = camera_to_rgb.map(|x| 1. / x.iter().sum::<f64>());
+	for (index, row) in camera_to_rgb.iter_mut().enumerate() {
+		*row = row.map(|x| x * white_balance_multiplier[index]);
+	}
+	let rgb_to_camera = transpose(pseudoinverse(camera_to_rgb));
+
+	raw_image.white_balance_multiplier = Some(white_balance_multiplier);
+	raw_image.camera_to_rgb = Some(camera_to_rgb);
+	raw_image.rgb_to_camera = Some(rgb_to_camera);
+
+	raw_image
+}
+
+#[allow(clippy::needless_range_loop)]
+fn pseudoinverse<const N: usize>(matrix: [[f64; 3]; N]) -> [[f64; 3]; N] {
+	let mut output_matrix = [[0.; 3]; N];
+	let mut work = [[0.; 6]; 3];
+
+	for i in 0..3 {
+		for j in 0..6 {
+			work[i][j] = if j == i + 3 { 1. } else { 0. };
+		}
+		for j in 0..3 {
+			for k in 0..N {
+				work[i][j] += matrix[k][i] * matrix[k][j];
+			}
+		}
+	}
+
+	for i in 0..3 {
+		let num = work[i][i];
+		for j in 0..6 {
+			work[i][j] /= num;
+		}
+		for k in 0..3 {
+			if k == i {
+				continue;
+			}
+			let num = work[k][i];
+			for j in 0..6 {
+				work[k][j] -= work[i][j] * num;
+			}
+		}
+	}
+
+	for i in 0..N {
+		for j in 0..3 {
+			output_matrix[i][j] = 0.;
+			for k in 0..3 {
+				output_matrix[i][j] += work[j][k + 3] * matrix[i][k];
+			}
+		}
+	}
+
+	output_matrix
+}
+
+fn transpose<const N: usize>(matrix: [[f64; 3]; N]) -> [[f64; N]; 3] {
+	let mut output_matrix = [[0.; N]; 3];
+
+	for (i, row) in matrix.iter().enumerate() {
+		for (j, &value) in row.iter().enumerate() {
+			output_matrix[j][i] = value;
+		}
+	}
+
+	output_matrix
 }

libraries/raw-rs/src/preprocessing/scale_colors.rs

@@ -1,44 +1,33 @@
 use crate::RawImage;
 
-const XYZ_TO_RGB: [[f64; 3]; 3] = [[0.412453, 0.357580, 0.180423], [0.212671, 0.715160, 0.072169], [0.019334, 0.119193, 0.950227]];
-
 pub fn scale_colors(mut raw_image: RawImage) -> RawImage {
-	if let Some(camera_to_xyz) = raw_image.camera_to_xyz {
-		let mut camera_to_rgb = [[0.; 3]; 3];
-		for i in 0..3 {
-			for j in 0..3 {
-				for k in 0..3 {
-					camera_to_rgb[i][j] += camera_to_xyz[i * 3 + k] * XYZ_TO_RGB[k][j];
-				}
-			}
-		}
+	let Some(mut white_balance_multiplier) = raw_image.white_balance_multiplier else {
+		return raw_image;
+	};
 
-		let mut white_balance_multiplier = camera_to_rgb.map(|x| 1. / x.iter().sum::<f64>());
+	if white_balance_multiplier[1] == 0. {
+		white_balance_multiplier[1] = 1.;
+	}
 
-		if white_balance_multiplier[1] == 0. {
-			white_balance_multiplier[1] = 1.;
-		}
+	// TODO: Move this at its correct location when highlights are implemented correctly.
+	let highlight = 0;
 
-		// TODO: Move this at its correct location when highlights are implemented correctly.
-		let highlight = 0;
+	let normalize_white_balance = if highlight == 0 {
+		white_balance_multiplier.iter().copied().fold(f64::INFINITY, f64::min)
+	} else {
+		white_balance_multiplier.iter().copied().fold(f64::NEG_INFINITY, f64::max)
+	};
 
-		let normalize_white_balance = if highlight == 0 {
-			white_balance_multiplier.iter().fold(f64::INFINITY, |a, &b| a.min(b))
-		} else {
-			white_balance_multiplier.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b))
-		};
+	let final_multiplier = if normalize_white_balance > 0.00001 && raw_image.maximum > 0 {
+		let scale_to_16bit_multiplier = u16::MAX as f64 / raw_image.maximum as f64;
+		white_balance_multiplier.map(|x| x / normalize_white_balance * scale_to_16bit_multiplier)
+	} else {
+		[1., 1., 1.]
+	};
 
-		let final_multiplier = if normalize_white_balance > 0.00001 && raw_image.maximum > 0 {
-			let scale_to_16bit_multiplier = u16::MAX as f64 / raw_image.maximum as f64;
-			white_balance_multiplier.map(|x| x / normalize_white_balance * scale_to_16bit_multiplier)
-		} else {
-			[1., 1., 1.]
-		};
-
-		for i in 0..(raw_image.height * raw_image.width) {
-			for (c, multiplier) in final_multiplier.iter().enumerate() {
-				raw_image.data[3 * i + c] = ((raw_image.data[3 * i + c] as f64) * multiplier).min(u16::MAX as f64).max(0.) as u16;
-			}
+	for i in 0..(raw_image.height * raw_image.width) {
+		for (c, multiplier) in final_multiplier.iter().enumerate() {
+			raw_image.data[3 * i + c] = ((raw_image.data[3 * i + c] as f64) * multiplier).min(u16::MAX as f64).max(0.) as u16;
 		}
 	}
 

libraries/raw-rs/tests/tests.rs

@@ -6,7 +6,6 @@ use raw_rs::RawImage;
 use image::codecs::png::{CompressionType, FilterType, PngEncoder};
 use image::{ColorType, ImageEncoder};
 use libraw::Processor;
-use palette::{LinSrgb, Srgb};
 use std::collections::HashMap;
 use std::fmt::Write;
 use std::fs::{read_dir, File};
@@ -69,16 +68,6 @@ fn test_images_match_with_libraw() {
 }
 
 fn store_image(path: &Path, suffix: &str, data: &mut [u8], width: usize, height: usize) {
-	if suffix == "raw_rs" {
-		for pixel in data.chunks_mut(3) {
-			let lin_srgb: LinSrgb<f64> = LinSrgb::new(pixel[0], pixel[1], pixel[2]).into_format();
-			let output: Srgb<u8> = Srgb::from_linear(lin_srgb);
-			pixel[0] = output.red;
-			pixel[1] = output.green;
-			pixel[2] = output.blue;
-		}
-	}
-
 	let mut output_path = PathBuf::new();
 	if let Some(parent) = path.parent() {
 		output_path.push(parent);