use std::collections::HashMap;
fn read(data: &[u8], index: usize, depth: u8) -> f32 {
match depth {
1 => data[index] as f32,
2 => u16::from_ne_bytes([data[index * 2], data[index * 2 + 1]]) as f32,
_ => f32::from_ne_bytes(data[index * 4..index * 4 + 4].try_into().unwrap()),
}
}
fn write(output: &mut [u8], index: usize, value: f32, depth: u8) {
match depth {
1 => output[index] = value.round_ties_even().clamp(0., 255.) as u8,
2 => output[index * 2..index * 2 + 2]
.copy_from_slice(&(value.round_ties_even().clamp(0., 65535.) as u16).to_ne_bytes()),
_ => output[index * 4..index * 4 + 4].copy_from_slice(&value.to_ne_bytes()),
}
}
#[allow(clippy::too_many_arguments)]
pub fn resize(
data: &[u8],
width: usize,
height: usize,
channels: usize,
depth: u8,
target_w: usize,
target_h: usize,
mode: u8,
xindices: &[Vec<usize>],
xweights: &[Vec<f32>],
yindices: &[Vec<usize>],
yweights: &[Vec<f32>],
) -> Result<Vec<u8>, String> {
if data.len()
!= width
.checked_mul(height)
.and_then(|v| v.checked_mul(channels))
.and_then(|v| v.checked_mul(depth as usize))
.ok_or("image size overflow")?
|| xindices.len()
!= if mode == 4 {
target_w * channels
} else {
target_w
}
|| yindices.len() != target_h
|| xweights.len() != xindices.len()
|| yweights.len() != target_h
|| xindices.iter().flatten().any(|&x| x >= width)
|| yindices.iter().flatten().any(|&y| y >= height)
|| xindices
.iter()
.zip(xweights)
.any(|(i, w)| i.len() != w.len())
|| yindices
.iter()
.zip(yweights)
.any(|(i, w)| i.len() != w.len())
{
return Err("invalid resize buffer or coefficient table".into());
}
let output_size = target_w
.checked_mul(target_h)
.and_then(|v| v.checked_mul(channels))
.and_then(|v| v.checked_mul(depth as usize))
.ok_or("output size overflow")?;
let mut output = vec![0; output_size];
if depth == 1 && [4, 5].contains(&mode) {
let bits = if mode == 4 { 8 } else { 7 };
let rounding = 1i32 << (bits - 1);
for y in 0..target_h {
for x in 0..target_w {
for c in 0..channels {
let coefficient = if mode == 4 { x * channels + c } else { x };
let vertical = |column: usize| -> i32 {
let a = data[(yindices[y][0] * width + xindices[coefficient][column])
* channels
+ c] as i32;
let b = data[(yindices[y][1] * width + xindices[coefficient][column])
* channels
+ c] as i32;
(a * yweights[y][0] as i32 + b * yweights[y][1] as i32 + rounding) >> bits
};
output[(y * target_w + x) * channels + c] =
((vertical(0) * xweights[coefficient][0] as i32
+ vertical(1) * xweights[coefficient][1] as i32
+ rounding)
>> bits) as u8;
}
}
}
return Ok(output);
}
let mut rows: HashMap<usize, Vec<f32>> = HashMap::new();
for y in 0..target_h {
rows.retain(|key, _| yindices[y].contains(key));
for &source_y in &yindices[y] {
rows.entry(source_y).or_insert_with(|| {
let mut row = vec![0f32; target_w * channels];
for x in 0..target_w {
for channel in 0..channels {
let mut sum = 0f32;
for (&index, &weight) in xindices[x].iter().zip(&xweights[x]) {
let pixel =
read(data, (source_y * width + index) * channels + channel, depth);
sum += pixel * weight;
}
row[x * channels + channel] = sum;
}
}
row
});
}
let vertical_rows: Vec<&[f32]> = yindices[y]
.iter()
.map(|index| rows[index].as_slice())
.collect();
for x in 0..target_w * channels {
if depth == 1 && mode == 1 {
let mut sum = 0i32;
for (tap, &weight) in yweights[y].iter().enumerate() {
sum += (((vertical_rows[tap][x] as i64 >> 4) * weight as i64) >> 16) as i32;
}
output[(y * target_w * channels) + x] = ((sum + 2) >> 2).clamp(0, 255) as u8;
} else if depth == 1 && mode == 6 && x < target_w * channels / 8 * 8 {
let mut sum = vertical_rows[3][x] * (yweights[y][3] / (1u32 << 22) as f32);
for index in (0..3).rev() {
sum = vertical_rows[index][x]
.mul_add(yweights[y][index] / (1u32 << 22) as f32, sum);
}
write(&mut output, y * target_w * channels + x, sum, depth);
} else if depth == 1 && [2, 6].contains(&mode) {
let mut sum = 0i64;
for (tap, &weight) in yweights[y].iter().enumerate() {
sum += vertical_rows[tap][x] as i64 * weight as i64;
}
output[y * target_w * channels + x] = ((sum + (1 << 21)) >> 22).clamp(0, 255) as u8;
} else {
let mut sum = 0f32;
for (tap, &weight) in yweights[y].iter().enumerate() {
sum += vertical_rows[tap][x] * weight;
}
if mode == 7 {
sum = (sum + 0.5).floor();
}
write(&mut output, y * target_w * channels + x, sum, depth);
}
}
}
Ok(output)
}
fn find(parents: &mut [usize], mut index: usize) -> usize {
while parents[index] != index {
parents[index] = parents[parents[index]];
index = parents[index];
}
index
}
pub fn components(
data: &[u8],
width: usize,
height: usize,
) -> Result<(Vec<u8>, Vec<[i32; 5]>), String> {
if width.checked_mul(height) != Some(data.len())
|| width > i32::MAX as usize
|| height > i32::MAX as usize
{
return Err("invalid component mask".into());
}
let mut labels = vec![0usize; data.len()];
let mut parents = vec![0usize];
for y in 0..height {
for x in 0..width {
if data[y * width + x] == 0 {
continue;
}
let mut roots = [0usize; 4];
let neighbors = [
if x > 0 { labels[y * width + x - 1] } else { 0 },
if y > 0 && x > 0 {
labels[(y - 1) * width + x - 1]
} else {
0
},
if y > 0 {
labels[(y - 1) * width + x]
} else {
0
},
if y > 0 && x + 1 < width {
labels[(y - 1) * width + x + 1]
} else {
0
},
];
let mut minimum = usize::MAX;
for (i, neighbor) in neighbors.into_iter().enumerate() {
if neighbor != 0 {
roots[i] = find(&mut parents, neighbor);
minimum = minimum.min(roots[i]);
}
}
let label = if minimum != usize::MAX {
for root in roots {
if root != 0 {
parents[root] = minimum;
}
}
minimum
} else {
let index = parents.len();
parents.push(index);
index
};
labels[y * width + x] = label;
}
}
let roots: Vec<usize> = (0..parents.len()).map(|v| find(&mut parents, v)).collect();
let mut mapping = vec![0usize; parents.len()];
let mut count = 1;
for y in (0..height).step_by(2) {
for x in (0..width).step_by(2) {
for dy in 0..2 {
for dx in 0..2 {
if y + dy < height && x + dx < width {
let root = roots[labels[(y + dy) * width + x + dx]];
if root != 0 && mapping[root] == 0 {
mapping[root] = count;
count += 1;
}
}
}
}
}
}
let mut stats = vec![[i32::MAX, i32::MAX, -1, -1, 0]; count];
let mut output = Vec::with_capacity(data.len() * 4);
for (index, label) in labels.into_iter().enumerate() {
let id = mapping[roots[label]];
output.extend_from_slice(&(id as i32).to_ne_bytes());
let x = (index % width) as i32;
let y = (index / width) as i32;
let stat = &mut stats[id];
stat[0] = stat[0].min(x);
stat[1] = stat[1].min(y);
stat[2] = stat[2].max(x);
stat[3] = stat[3].max(y);
stat[4] += 1;
}
for stat in &mut stats {
if stat[4] == 0 {
*stat = [-1, -1, 0, 0, 0];
} else {
stat[2] -= stat[0] - 1;
stat[3] -= stat[1] - 1;
}
}
Ok((output, stats))
}
pub fn contours(
data: &[u8],
width: usize,
height: usize,
external: bool,
) -> Result<Vec<Vec<[i32; 2]>>, String> {
if width.checked_mul(height) != Some(data.len()) {
return Err("invalid contour mask".into());
}
let stride = width + 2;
let mut pixels = vec![0i16; stride * (height + 2)];
for y in 0..height {
for x in 0..width {
pixels[(y + 1) * stride + x + 1] = i16::from(data[y * width + x] != 0);
}
}
let deltas: [isize; 8] = [
1,
1 - stride as isize,
-(stride as isize),
-1 - stride as isize,
-1,
stride as isize - 1,
stride as isize,
stride as isize + 1,
];
let mut result = Vec::new();
for y in 1..=height {
let mut previous = 0i16;
let mut last_boundary = y * stride;
for x in 1..=width {
let index = y * stride + x;
let current = pixels[index];
if current == previous {
continue;
}
let outer = previous == 0 && current == 1;
let hole = current == 0 && previous >= 1;
if !outer && !hole {
previous = current;
if previous & -2 != 0 {
last_boundary = index;
}
continue;
}
if hole && previous & -2 != 0 {
last_boundary = index - 1;
}
if external && (hole || pixels[last_boundary] > 0) {
previous = current;
continue;
}
let start = index - usize::from(hole);
let mut direction = if hole { 0 } else { 4 };
let end = direction;
let mut first;
loop {
direction = (direction + 7) & 7;
first = (start as isize + deltas[direction]) as usize;
if pixels[first] != 0 || direction == end {
break;
}
}
let mut points = Vec::new();
if direction == end {
pixels[start] = -126;
points.push([(start % stride) as i32 - 1, (start / stride) as i32 - 1]);
} else {
let mut position = start;
let mut previous_direction = direction ^ 4;
let mut terminated = false;
for _ in 0..pixels.len() * 8 {
let end = direction;
let mut next_direction = direction + 1;
let mut following;
loop {
following = (position as isize + deltas[next_direction & 7]) as usize;
if pixels[following] != 0 {
break;
}
next_direction += 1;
}
direction = next_direction & 7;
if (direction as u32).wrapping_sub(1) < end as u32 {
pixels[position] = -126;
} else if pixels[position] == 1 {
pixels[position] = 2;
}
if direction != previous_direction {
points.push([
(position % stride) as i32 - 1,
(position / stride) as i32 - 1,
]);
previous_direction = direction;
}
if following == start && position == first {
terminated = true;
break;
}
position = following;
direction = (direction + 4) & 7;
}
if !terminated {
return Err("contour traversal did not terminate".into());
}
}
result.push(points);
last_boundary = start;
previous = pixels[index];
}
}
result.reverse();
Ok(result)
}
fn half_round(value: f64) -> f32 {
let bits = value.to_bits();
let exponent = ((bits >> 52) & 2047) as i32 - 1023;
if exponent < -14 {
return ((value * 16777216.).round_ties_even() / 16777216.) as f32;
}
if exponent >= 16 {
return if value.is_sign_negative() {
f32::NEG_INFINITY
} else {
f32::INFINITY
};
}
let mask = (1u64 << 42) - 1;
let remainder = bits & mask;
let mut rounded = bits & !mask;
if remainder > 1u64 << 41 || (remainder == 1u64 << 41 && (rounded >> 42) & 1 != 0) {
rounded += 1u64 << 42;
}
f64::from_bits(rounded) as f32
}
fn half_fma(a: f32, b: f32, c: f32) -> f32 {
half_round(a as f64 * b as f64 + c as f64)
}
fn cubic_weights(x: f32, half: bool) -> [f32; 4] {
if half {
let x = half_round(x as f64);
let a2 = half_round((x * x) as f64);
let b = half_round((1. - x) as f64);
let b2 = half_round((b * b) as f64);
let w0 = half_round((-0.75 * half_round((x * b2) as f64)) as f64);
let w3 = half_round((-0.75 * half_round((a2 * b) as f64)) as f64);
let w1 = half_fma(a2, half_fma(1.25, x, -2.25), 1.);
let w2 = half_round((half_round((half_round((1. - w0) as f64) - w1) as f64) - w3) as f64);
[w0, w1, w2, w3]
} else {
let a2 = x * x;
let b = 1. - x;
let w0 = -0.75 * (x * (b * b));
let w3 = -0.75 * (a2 * b);
let w1 = a2.mul_add(1.25f32.mul_add(x, -2.25), 1.);
[w0, w1, 1. - w0 - w1 - w3, w3]
}
}
#[allow(clippy::too_many_arguments)]
fn sample(
data: &[u8],
width: usize,
height: usize,
channels: usize,
depth: u8,
x: i64,
y: i64,
c: usize,
replicate: bool,
value: f32,
) -> f32 {
if !replicate && (x < 0 || y < 0 || x >= width as i64 || y >= height as i64) {
return value;
}
let x = x.clamp(0, width as i64 - 1) as usize;
let y = y.clamp(0, height as i64 - 1) as usize;
read(data, (y * width + x) * channels + c, depth)
}
#[allow(clippy::too_many_arguments)]
pub fn warp(
data: &[u8],
width: usize,
height: usize,
channels: usize,
depth: u8,
target_w: usize,
target_h: usize,
mode: u8,
replicate: bool,
value: f32,
inverse: &[f64; 9],
affine: bool,
) -> Result<Vec<u8>, String> {
if width == 0
|| height == 0
|| data.len()
!= width
.checked_mul(height)
.and_then(|v| v.checked_mul(channels))
.and_then(|v| v.checked_mul(depth as usize))
.ok_or("image size overflow")?
{
return Err("invalid warp image".into());
}
let length = target_w
.checked_mul(target_h)
.and_then(|v| v.checked_mul(channels))
.and_then(|v| v.checked_mul(depth as usize))
.ok_or("output size overflow")?;
let mut output = vec![0u8; length];
let half = depth == 1 && channels != 2;
let m: Vec<f32> = inverse.iter().map(|&v| v as f32).collect();
for y in 0..target_h {
let row: Vec<f32> = (0..3)
.map(|i| {
if depth == 4 {
(y as f32).mul_add(m[i * 3 + 1], m[i * 3 + 2])
} else {
(y as f64 * inverse[i * 3 + 1] + inverse[i * 3 + 2]) as f32
}
})
.collect();
let cubic_row: Vec<f32> = (0..3)
.map(|i| (y as f32).mul_add(m[i * 3 + 1], m[i * 3 + 2]))
.collect();
for x in 0..target_w {
let (sx, sy) = if mode == 2 {
if affine {
(
(x as f32).mul_add(m[0], cubic_row[0]),
(x as f32).mul_add(m[3], cubic_row[1]),
)
} else {
let d = cubic_row[2] as f64 + x as f64 * m[6] as f64;
if d == 0. {
(0., 0.)
} else {
(
((cubic_row[0] as f64 + x as f64 * m[0] as f64) / d) as f32,
((cubic_row[1] as f64 + x as f64 * m[3] as f64) / d) as f32,
)
}
}
} else if x < target_w / 8 * 8 {
let d = (x as f32).mul_add(m[6], row[2]);
if d == 0. {
(0., 0.)
} else {
(
(x as f32).mul_add(m[0], row[0]) / d,
(x as f32).mul_add(m[3], row[1]) / d,
)
}
} else {
let d = (x as f64 * inverse[6] + y as f64 * inverse[7] + inverse[8]) as f32;
if d == 0. {
(0., 0.)
} else {
(
((x as f64 * inverse[0] + y as f64 * inverse[1] + inverse[2]) / d as f64)
as f32,
((x as f64 * inverse[3] + y as f64 * inverse[4] + inverse[5]) / d as f64)
as f32,
)
}
};
let bx = sx.floor() as i64;
let by = sy.floor() as i64;
let alpha = sx - bx as f32;
let beta = sy - by as f32;
let wx = cubic_weights(alpha, half);
let wy = cubic_weights(beta, half);
for c in 0..channels {
let pixel = |ix: i64, iy: i64| {
sample(
data, width, height, channels, depth, ix, iy, c, replicate, value,
)
};
let result = if mode == 0 {
pixel(sx.round_ties_even() as i64, sy.round_ties_even() as i64)
} else if mode == 1 {
let p00 = pixel(bx, by);
let p01 = pixel(bx + 1, by);
let p10 = pixel(bx, by + 1);
let p11 = pixel(bx + 1, by + 1);
let top = alpha.mul_add(p01 - p00, p00);
let bottom = alpha.mul_add(p11 - p10, p10);
beta.mul_add(bottom - top, top)
} else {
let mut sum = 0f32;
for (j, &yweight) in wy.iter().enumerate() {
let mut horizontal = pixel(bx - 1, by + j as i64 - 1) * wx[0];
if half {
horizontal = half_round(horizontal as f64);
}
for (i, &xweight) in wx.iter().enumerate().skip(1) {
let p = pixel(bx + i as i64 - 1, by + j as i64 - 1);
horizontal = if half {
half_fma(p, xweight, horizontal)
} else {
p.mul_add(xweight, horizontal)
};
}
sum = if half {
half_fma(horizontal, yweight, sum)
} else {
horizontal.mul_add(yweight, sum)
};
}
sum
};
write(
&mut output,
(y * target_w + x) * channels + c,
result,
depth,
);
}
}
}
Ok(output)
}
pub fn minimum_rectangle(hull: &[[f32; 2]]) -> [f32; 5] {
let n = hull.len();
if n == 0 {
return [0., 0., 0., 0., -90.];
}
if n == 1 {
return [hull[0][0], hull[0][1], 0., 0., -90.];
}
if n == 2 {
let dx = (hull[0][0] - hull[1][0]) as f64;
let dy = (hull[0][1] - hull[1][1]) as f64;
let mut w = 0.;
let mut h = dx.hypot(dy) as f32;
let mut angle = -90.;
if dx == 0. {
std::mem::swap(&mut w, &mut h);
} else if dy < 0. {
angle = dy.atan2(dx).to_degrees() as f32;
std::mem::swap(&mut w, &mut h);
} else if dy > 0. {
angle = -dx.atan2(dy).to_degrees() as f32;
}
return [
(hull[0][0] + hull[1][0]) * 0.5,
(hull[0][1] + hull[1][1]) * 0.5,
w,
h,
angle,
];
}
let vectors: Vec<[f32; 2]> = (0..n)
.map(|i| {
[
hull[(i + 1) % n][0] - hull[i][0],
hull[(i + 1) % n][1] - hull[i][1],
]
})
.collect();
let lengths: Vec<f32> = vectors
.iter()
.map(|v| (1. / (v[0] as f64).hypot(v[1] as f64)) as f32)
.collect();
let extremum = |axis: usize, maximum: bool| {
let mut idx = 0;
for i in 1..n {
if if maximum {
hull[i][axis] > hull[idx][axis]
} else {
hull[i][axis] < hull[idx][axis]
} {
idx = i;
}
}
idx
};
let mut seq = [
extremum(1, false),
extremum(0, true),
extremum(1, true),
extremum(0, false),
];
let mut area = f32::INFINITY;
let (mut left, mut a, mut b, mut width, mut height, mut bottom) =
(0, 0f32, 0f32, 0f32, 0f32, 0);
for _ in 0..n {
let v0 = vectors[seq[0]];
let v1 = vectors[seq[1]];
let v2 = vectors[seq[2]];
let v3 = vectors[seq[3]];
let rot = [v0, [v1[1], -v1[0]], [-v2[0], -v2[1]], [-v3[1], v3[0]]];
let mut main = 0;
for i in 1..4 {
if rot[i][1].mul_add(rot[main][0], -rot[i][0] * rot[main][1]) < 0. {
main = i;
}
}
let lead_x = vectors[seq[main]][0] * lengths[seq[main]];
let lead_y = vectors[seq[main]][1] * lengths[seq[main]];
let (ba, bb) = [
(lead_x, lead_y),
(lead_y, -lead_x),
(-lead_x, -lead_y),
(-lead_y, lead_x),
][main];
seq[main] = (seq[main] + 1) % n;
let dx = hull[seq[1]][0] - hull[seq[3]][0];
let dy = hull[seq[1]][1] - hull[seq[3]][1];
let bw = dx.mul_add(ba, dy * bb);
let dx = hull[seq[2]][0] - hull[seq[0]][0];
let dy = hull[seq[2]][1] - hull[seq[0]][1];
let bh = (-dx).mul_add(bb, dy * ba);
if bw * bh <= area {
area = bw * bh;
left = seq[3];
a = ba;
b = bb;
width = bw;
height = bh;
bottom = seq[0];
}
}
let c1 = a.mul_add(hull[left][0], hull[left][1] * b);
let c2 = (-b).mul_add(hull[bottom][0], hull[bottom][1] * a);
let inverse = 1. / a.mul_add(a, b * b);
let px = c1.mul_add(a, -c2 * b) * inverse;
let py = a.mul_add(c2, b * c1) * inverse;
let edge1 = [a * width, b * width];
let edge2 = [-b * height, a * height];
let mut w = (edge2[0] as f64).hypot(edge2[1] as f64) as f32;
let mut h = (edge1[0] as f64).hypot(edge1[1] as f64) as f32;
let angle = if edge1[0] == 0. && edge1[1] > 0. {
std::mem::swap(&mut w, &mut h);
-90.
} else {
-(edge1[0] as f64).atan2(edge1[1] as f64).to_degrees() as f32
};
[
px + (edge1[0] + edge2[0]) * 0.5,
py + (edge1[1] + edge2[1]) * 0.5,
w,
h,
angle,
]
}
pub fn morphology(
data: &[u8],
width: usize,
height: usize,
kw: usize,
kh: usize,
dilate: bool,
) -> Result<Vec<u8>, String> {
if width.checked_mul(height) != Some(data.len()) || kw == 0 || kh == 0 {
return Err("invalid morphology buffer".into());
}
let border = if dilate { 0 } else { 255 };
let pass = |input: &[u8], vertical: bool, kernel: usize| {
let length = if vertical { height } else { width };
let lines = if vertical { width } else { height };
let mut output = vec![0u8; data.len()];
for line in 0..lines {
let mut queue = std::collections::VecDeque::<(usize, u8)>::new();
for p in 0..length + kernel - 1 {
let coordinate = p as isize - (kernel / 2) as isize;
let value = if coordinate < 0 || coordinate >= length as isize {
border
} else {
input[if vertical {
coordinate as usize * width + line
} else {
line * width + coordinate as usize
}]
};
while queue
.back()
.is_some_and(|v| if dilate { v.1 <= value } else { v.1 >= value })
{
queue.pop_back();
}
queue.push_back((p, value));
while queue.front().is_some_and(|v| v.0 + kernel <= p) {
queue.pop_front();
}
if p + 1 >= kernel {
let position = p + 1 - kernel;
output[if vertical {
position * width + line
} else {
line * width + position
}] = queue.front().unwrap().1;
}
}
}
output
};
Ok(pass(&pass(data, false, kw), true, kh))
}
pub fn convex_hull(points: &[[f32; 2]]) -> Vec<[f32; 2]> {
let mut values: Vec<([f32; 2], usize)> = points
.iter()
.copied()
.enumerate()
.map(|(i, p)| (p, i))
.collect();
values.sort_unstable_by(|a, b| {
a.0[0]
.total_cmp(&b.0[0])
.then_with(|| a.0[1].total_cmp(&b.0[1]))
.then(a.1.cmp(&b.1))
});
values.dedup_by(|a, b| a.0 == b.0);
if values.len() < 2 {
return values.into_iter().map(|v| v.0).collect();
}
let cross = |a: ([f32; 2], usize), b: ([f32; 2], usize), c: ([f32; 2], usize)| {
(b.0[0] - a.0[0]) as f64 * (c.0[1] - a.0[1]) as f64
- (b.0[1] - a.0[1]) as f64 * (c.0[0] - a.0[0]) as f64
};
let mut lower = Vec::new();
let mut upper = Vec::new();
for &point in &values {
while lower.len() >= 2 && cross(lower[lower.len() - 2], lower[lower.len() - 1], point) <= 0.
{
lower.pop();
}
lower.push(point);
}
for &point in values.iter().rev() {
while upper.len() >= 2 && cross(upper[upper.len() - 2], upper[upper.len() - 1], point) <= 0.
{
upper.pop();
}
upper.push(point);
}
lower.pop();
upper.pop();
lower.extend(upper);
let n = lower.len();
let mut start = 0;
for i in 1..n {
if lower[i].0[0] > lower[start].0[0]
|| (lower[i].0[0] == lower[start].0[0] && lower[i].0[1] > lower[start].0[1])
{
start = i;
}
}
lower.rotate_left(start);
if n >= 3 {
for increasing in [true, false] {
let start = if increasing {
lower.iter().enumerate().min_by_key(|(_, v)| v.1).unwrap().0
} else {
lower.iter().enumerate().max_by_key(|(_, v)| v.1).unwrap().0
};
if (0..n - 1)
.all(|i| (lower[(start + i) % n].1 < lower[(start + i + 1) % n].1) == increasing)
{
lower.rotate_left(start);
break;
}
}
}
lower.into_iter().map(|v| v.0).collect()
}
pub fn linear_solve(matrix: &[Vec<f64>], rhs: &[f64]) -> Option<Vec<f64>> {
let n = matrix.len();
if n == 0 || n > 16 || rhs.len() != n || matrix.iter().any(|row| row.len() != n) {
return None;
}
let mut a = matrix.to_vec();
let mut b = rhs.to_vec();
for i in 0..n {
let mut pivot = i;
for j in i + 1..n {
if a[j][i].abs() > a[pivot][i].abs() {
pivot = j;
}
}
if a[pivot][i].abs() < f64::EPSILON * 100. {
return None;
}
if pivot != i {
a.swap(i, pivot);
b.swap(i, pivot);
}
let reciprocal = -1. / a[i][i];
let (prior, following) = a.split_at_mut(i + 1);
let pivot_row = &prior[i];
for (offset, row) in following.iter_mut().enumerate() {
let alpha = row[i] * reciprocal;
for (value, &coefficient) in row[i + 1..].iter_mut().zip(&pivot_row[i + 1..]) {
*value = alpha.mul_add(coefficient, *value);
}
let j = i + 1 + offset;
b[j] = alpha.mul_add(b[i], b[j]);
}
}
for i in (0..n).rev() {
let mut sum = b[i];
for (k, &value) in b.iter().enumerate().skip(i + 1) {
sum = (-a[i][k]).mul_add(value, sum);
}
b[i] = sum / a[i][i];
}
Some(b)
}