From 8c9dde4c210b81a3df3b767033f9c308527f71ad Mon Sep 17 00:00:00 2001 From: Martin Bruse Date: Mon, 16 Dec 2024 19:28:25 +0100 Subject: [PATCH] Added more bits and pieces to the spline rendering. --- jxl/src/error.rs | 10 +- jxl/src/features/spline.rs | 406 +++++++++++++++++++++++++++++++++++-- jxl/src/util/test.rs | 4 +- 3 files changed, 402 insertions(+), 18 deletions(-) diff --git a/jxl/src/error.rs b/jxl/src/error.rs index 2585959..641e6f2 100644 --- a/jxl/src/error.rs +++ b/jxl/src/error.rs @@ -7,7 +7,9 @@ use std::collections::TryReserveError; use thiserror::Error; -use crate::{entropy_coding::huffman::HUFFMAN_MAX_BITS, image::DataTypeTag}; +use crate::{ + entropy_coding::huffman::HUFFMAN_MAX_BITS, features::spline::Point, image::DataTypeTag, +}; #[derive(Error, Debug)] pub enum Error { @@ -111,18 +113,22 @@ pub enum Error { InvalidPredictor(u32), #[error("Invalid modular mode property: {0}")] InvalidProperty(u32), + #[error("Point list is empty")] + PointListEmpty, #[error("Too large area for spline: {0}, limit is {1}")] SplinesAreaTooLarge(u64, u64), #[error("Too large manhattan_distance reached: {0}, limit is {1}")] SplinesDistanceTooLarge(u64, u64), #[error("Too many splines: {0}, limit is {1}")] SplinesTooMany(u32, u32), + #[error("Spline has adjacent coinciding control points: point[{0}]: {1:?}, point[{2}]: {3:?}")] + SplineAdjacentCoincidingControlPoints(u32, Point, u32, Point), #[error("Too many control points for splines: {0}, limit is {1}")] SplinesTooManyControlPoints(u32, u32), #[error( "Spline point outside valid bounds: coordinates: {0:?}, out of bounds: {1}, bounds: {2:?}" )] - SplinesPointOutOfRange((i32, i32), i32, std::ops::Range), + SplinesPointOutOfRange(Point, i32, std::ops::Range), #[error("Spline coordinates out of bounds: {0}, limit is {1}")] SplinesCoordinatesLimit(i32, i32), #[error("Spline delta-delta is out of bounds: {0}, limit is {1}")] diff --git a/jxl/src/features/spline.rs b/jxl/src/features/spline.rs index 591d1b7..bae01e1 100644 --- a/jxl/src/features/spline.rs +++ b/jxl/src/features/spline.rs @@ -5,7 +5,10 @@ // TODO(firsching): remove once we use this! #![allow(dead_code)] -use std::f32::consts::FRAC_1_SQRT_2; +use std::{ + f32::consts::FRAC_1_SQRT_2, + ops::{self, Index}, +}; use crate::{ bit_reader::BitReader, @@ -26,7 +29,7 @@ const CONTROL_POINTS_CONTEXT: usize = 4; const DCT_CONTEXT: usize = 5; const NUM_SPLINE_CONTEXTS: usize = 6; -#[derive(Debug, Clone, Copy)] +#[derive(Debug, Clone, Copy, Default)] pub struct Point { pub x: f32, pub y: f32, @@ -36,6 +39,9 @@ impl Point { fn new(x: f32, y: f32) -> Self { Point { x, y } } + fn abs(&self) -> f32 { + self.x.hypot(self.y) + } } impl PartialEq for Point { @@ -44,6 +50,47 @@ impl PartialEq for Point { } } +impl ops::Add for Point { + type Output = Point; + fn add(self, rhs: Point) -> Point { + Point { + x: self.x + rhs.x, + y: self.y + rhs.y, + } + } +} + +impl ops::Sub for Point { + type Output = Point; + fn sub(self, rhs: Point) -> Point { + Point { + x: self.x - rhs.x, + y: self.y - rhs.y, + } + } +} + +impl ops::Mul for Point { + type Output = Point; + fn mul(self, rhs: f32) -> Point { + Point { + x: self.x * rhs, + y: self.y * rhs, + } + } +} + +impl ops::Div for Point { + type Output = Point; + fn div(self, rhs: f32) -> Point { + let inv = 1.0 / rhs; + Point { + x: self.x * inv, + y: self.y * inv, + } + } +} + pub struct Spline { control_points: Vec, // X, Y, B. @@ -55,6 +102,27 @@ pub struct Spline { estimated_area_reached: u64, } +impl Spline { + pub fn validate_adjacent_point_coincidence(&self) -> Result<()> { + if let Some(item) = self + .control_points + .iter() + .enumerate() + .find(|(index, point)| { + index + 1 < self.control_points.len() && self.control_points[index + 1] == **point + }) + { + return Err(Error::SplineAdjacentCoincidingControlPoints( + item.0 as u32, + *item.1, + (item.0 + 1) as u32, + self.control_points[item.0 + 1], + )); + } + Ok(()) + } +} + #[derive(Debug, Default, Clone)] pub struct QuantizedSpline { // Double delta-encoded. @@ -76,16 +144,34 @@ fn validate_spline_point_pos(x: T, y: T) -> Result<( let yi = y.to_i32().unwrap(); let ok_range = -(1i32 << 23)..(1i32 << 23); if !ok_range.contains(&xi) { - return Err(Error::SplinesPointOutOfRange((xi, yi), xi, ok_range)); + return Err(Error::SplinesPointOutOfRange( + Point { + x: xi as f32, + y: yi as f32, + }, + xi, + ok_range, + )); } if !ok_range.contains(&yi) { - return Err(Error::SplinesPointOutOfRange((xi, yi), yi, ok_range)); + return Err(Error::SplinesPointOutOfRange( + Point { + x: xi as f32, + y: yi as f32, + }, + yi, + ok_range, + )); } Ok(()) } const CHANNEL_WEIGHT: [f32; 4] = [0.0042, 0.075, 0.07, 0.3333]; +fn area_limit(image_size: u64) -> u64 { + (1024 * image_size + (1u64 << 32)).min(1u64 << 42) +} + impl QuantizedSpline { #[instrument(level = "debug", skip(br), ret, err)] pub fn read( @@ -144,7 +230,7 @@ impl QuantizedSpline { y_to_b: f32, image_size: u64, ) -> Result { - let area_limit = (1024 * image_size + (1u64 << 32)).min(1u64 << 42); + let area_limit = area_limit(image_size); let mut result = Spline { control_points: Vec::with_capacity(self.control_points.len() + 1), @@ -237,15 +323,6 @@ impl QuantizedSpline { result.estimated_area_reached = width_estimate * manhattan_distance; - // TODO: move this check to the outside, using the returned estimated_area_reached, - // making the caller keep track of the total estimated_area_readed - //if result.total_estimated_area_reached > area_limit { - // return Err(Error::SplinesAreaTooLarge( - // *total_estimated_area_reached, - // area_limit, - // )); - //} - Ok(result) } } @@ -270,11 +347,213 @@ pub struct Splines { segment_y_start: Vec, } +struct ExtendedPoints<'a> { + first_element: Point, + last_element: Point, + backend: &'a Vec, +} + +impl ExtendedPoints<'_> { + fn len(&self) -> usize { + self.backend.len() + 2 + } + fn iter(&self) -> ExtendedPointsIterator { + ExtendedPointsIterator { + index: 0, + backend: self, + } + } +} + +impl Index for ExtendedPoints<'_> { + type Output = Point; + fn index(&self, index: usize) -> &Self::Output { + if index == 0 { + &self.first_element + } else if index == self.len() - 1 { + &self.last_element + } else { + &self.backend[index - 1] + } + } +} + +struct ExtendedPointsIterator<'a> { + index: usize, + backend: &'a ExtendedPoints<'a>, +} + +impl Iterator for ExtendedPointsIterator<'_> { + type Item = Point; + + fn next(&mut self) -> Option { + if self.index < self.backend.len() { + let result = Some(self.backend[self.index]); + self.index += 1; + result + } else { + None + } + } +} + +fn draw_centripetal_catmull_rom_spline(points: &Vec) -> Vec { + if points.is_empty() { + return [].to_vec(); + } + if points.len() == 1 { + return vec![points[0]]; + } + const NUM_POINTS: usize = 16; + // Create a view of points with one prepended and one appended point. + let extended_points = ExtendedPoints { + first_element: points[0] + (points[0] - points[1]), + last_element: points[points.len() - 1] + + (points[points.len() - 1] - points[points.len() - 2]), + backend: points, + }; + // Precompute the deltas used in the loop. + let deltas: Vec = extended_points + .iter() + .enumerate() + .map(|(index, point)| { + if index + 1 < extended_points.len() { + (extended_points[index + 1] - point).abs().sqrt() + } else { + 0.0 + } + }) + .collect(); + let mut result = Vec::::with_capacity((points.len() - 1) * NUM_POINTS + 1); + for start in 0..(extended_points.len() - 3) { + result.push(extended_points[start + 1]); + let mut t = [0.0; 4]; + for k in 0..3 { + // TODO(from libjxl): Restrict d[k] with reasonable limit and spec it. + t[k + 1] = t[k] + deltas[start + k]; + } + for i in 1..NUM_POINTS { + let tt = deltas[start] + ((i as f32) / (NUM_POINTS as f32)) * deltas[start + 1]; + let mut a = [Point::default(); 3]; + for k in 0..3 { + // TODO(from libjxl): Reciprocal multiplication would be faster. + a[k] = extended_points[start + k] + + (extended_points[start + k + 1] - extended_points[start + k]) + * ((tt - t[k]) / deltas[start + k]); + } + let mut b = [Point::default(); 2]; + for k in 0..2 { + b[k] = a[k] + + (a[k + 1] - a[k]) + * ((tt - t[k]) / (deltas[start + k] + deltas[start + k + 1])); + } + result.push(b[0] + (b[1] - b[0]) * ((tt - t[1]) / deltas[start + 1])); + } + } + result.push(extended_points[extended_points.len() - 2]); + result +} + +fn for_each_equally_spaced_point( + points: &[Point], + desired_distance: f32, + mut f: F, +) { + if points.is_empty() { + return; + } + let mut accumulated_distance = 0.0; + f(points[0], desired_distance); + if points.len() == 1 { + return; + } + for index in 0..(points.len() - 1) { + let mut current = points[index]; + let next = points[index + 1]; + let segment = next - current; + let segment_length = segment.abs(); + let unit_step = segment / segment_length; + if accumulated_distance + segment_length >= desired_distance { + current = current + unit_step * (desired_distance - accumulated_distance); + f(current, desired_distance); + accumulated_distance -= desired_distance; + } + accumulated_distance += segment_length; + while accumulated_distance >= desired_distance { + current = current + unit_step * desired_distance; + f(current, desired_distance); + accumulated_distance -= desired_distance; + } + } + f(points[points.len() - 1], accumulated_distance); +} + +const DESIRED_RENDERING_DISTANCE: f32 = 1.0; + impl Splines { fn has_any(&self) -> bool { !self.splines.is_empty() } + // TODO(zond): Add color correlation as parameter. + pub fn initialize_draw_cache(&mut self, image_xsize: u64, image_ysize: u64) -> Result<()> { + self.segments.clear(); + self.segment_indices.clear(); + self.segment_y_start.clear(); + // let mut segments_by_y = Vec::new(); + // let mut intermediate_points = Vec::new(); + let mut total_estimated_area_reached = 0u64; + let mut splines = Vec::new(); + // TODO(zond): Use color correlation here. + let y_to_x = 0.0; + let y_to_b = 0.0; + let area_limit = area_limit(image_xsize * image_ysize); + for (index, qspline) in self.splines.iter().enumerate() { + let spline = qspline.dequantize( + &self.starting_points[index], + self.quantization_adjustment, + y_to_x, + y_to_b, + image_xsize * image_ysize, + )?; + total_estimated_area_reached += spline.estimated_area_reached; + if total_estimated_area_reached > area_limit { + return Err(Error::SplinesAreaTooLarge( + total_estimated_area_reached, + area_limit, + )); + } + spline.validate_adjacent_point_coincidence()?; + splines.push(spline); + } + + if total_estimated_area_reached + > (8 * image_xsize * image_ysize + (1u64 << 25)).min(1u64 << 30) + { + warn!( + "Large total_estimated_area_reached, expect slower decoding:{}", + total_estimated_area_reached + ); + } + + for spline in splines { + let mut points_to_draw = Vec::<(Point, f32)>::new(); + let intermediate_points = draw_centripetal_catmull_rom_spline(&spline.control_points); + for_each_equally_spaced_point( + &intermediate_points, + DESIRED_RENDERING_DISTANCE, + |p, d| points_to_draw.push((p, d)), + ); + let length = (points_to_draw.len() - 2) as f32 * DESIRED_RENDERING_DISTANCE + + points_to_draw[points_to_draw.len() - 1].1; + if length <= 0.0 { + continue; + } + } + + todo!("finish translating this function from C++"); + } + #[instrument(level = "debug", skip(br), ret, err)] pub fn read(br: &mut BitReader, num_pixels: u32) -> Result { trace!(pos = br.total_bits_read()); @@ -347,7 +626,10 @@ impl Splines { mod test_splines { use crate::{error::Error, util::test::assert_all_almost_eq}; - use super::{Point, QuantizedSpline, Spline}; + use super::{ + draw_centripetal_catmull_rom_spline, for_each_equally_spaced_point, Point, QuantizedSpline, + Spline, + }; #[test] fn dequantize() -> Result<(), Error> { @@ -626,4 +908,98 @@ mod test_splines { } Ok(()) } + + #[test] + fn draw_spline() -> Result<(), Error> { + let control_points = vec![Point { x: 1.0, y: 2.0 }, Point { x: 4.0, y: 3.0 }]; + let want_result = [ + Point { x: 1.0, y: 2.0 }, + Point { + x: 1.1875, + y: 2.0625, + }, + Point { x: 1.375, y: 2.125 }, + Point { + x: 1.5625, + y: 2.1875, + }, + Point { x: 1.75, y: 2.25 }, + Point { + x: 1.9375, + y: 2.3125, + }, + Point { x: 2.125, y: 2.375 }, + Point { + x: 2.3125, + y: 2.4375, + }, + Point { x: 2.5, y: 2.5 }, + Point { + x: 2.6875, + y: 2.5625, + }, + Point { x: 2.875, y: 2.625 }, + Point { + x: 3.0625, + y: 2.6875, + }, + Point { x: 3.25, y: 2.75 }, + Point { + x: 3.4375, + y: 2.8125, + }, + Point { x: 3.625, y: 2.875 }, + Point { + x: 3.8125, + y: 2.9375, + }, + Point { x: 4.0, y: 3.0 }, + ]; + let got_result = draw_centripetal_catmull_rom_spline(&control_points); + assert_all_almost_eq!( + got_result.iter().map(|p| p.x).collect::>(), + want_result.iter().map(|p| p.x).collect::>(), + 1e-6, + ); + Ok(()) + } + + #[test] + fn equally_spaced_points() -> Result<(), Error> { + let desired_rendering_distance = 10.0f32; + let segments = [ + Point { x: 0.0, y: 0.0 }, + Point { x: 5.0, y: 0.0 }, + Point { x: 35.0, y: 0.0 }, + Point { x: 35.0, y: 10.0 }, + ]; + let want_results = [ + (Point { x: 0.0, y: 0.0 }, desired_rendering_distance), + (Point { x: 10.0, y: 0.0 }, desired_rendering_distance), + (Point { x: 20.0, y: 0.0 }, desired_rendering_distance), + (Point { x: 30.0, y: 0.0 }, desired_rendering_distance), + (Point { x: 35.0, y: 5.0 }, desired_rendering_distance), + (Point { x: 35.0, y: 10.0 }, 5.0f32), + ]; + let mut got_results = Vec::<(Point, f32)>::new(); + for_each_equally_spaced_point(&segments, desired_rendering_distance, |p, d| { + got_results.push((p, d)) + }); + assert_all_almost_eq!( + got_results.iter().map(|(p, _)| p.x).collect::>(), + want_results.iter().map(|(p, _)| p.x).collect::>(), + 1e-9 + ); + assert_all_almost_eq!( + got_results.iter().map(|(p, _)| p.y).collect::>(), + want_results.iter().map(|(p, _)| p.y).collect::>(), + 1e-9 + ); + assert_all_almost_eq!( + got_results.iter().map(|(_, d)| *d).collect::>(), + want_results.iter().map(|(_, d)| *d).collect::>(), + 1e-9 + ); + Ok(()) + } } diff --git a/jxl/src/util/test.rs b/jxl/src/util/test.rs index a6e8d13..1068435 100644 --- a/jxl/src/util/test.rs +++ b/jxl/src/util/test.rs @@ -28,7 +28,9 @@ pub(crate) use assert_almost_eq; macro_rules! assert_all_almost_eq { ($left:expr, $right:expr, $max_error:expr $(,)?) => { let (left_val, right_val, max_error) = (&$left, &$right, &$max_error); - assert_eq!(left_val.len(), right_val.len()); + if left_val.len() != right_val.len() { + panic!("assertion failed: `(left ≈ right)`\n left.len(): `{}`,\n right.len(): `{}`", left_val.len(), right_val.len()); + } for index in 0..left_val.len() { match $crate::util::test::abs_delta(left_val[index], right_val[index]).partial_cmp(max_error) { Some(std::cmp::Ordering::Greater) | None => panic!(