mirror of
https://github.com/storytold/lightcraft.git
synced 2026-10-14 19:27:01 -04:00
lensdb: lens distortion (poly3/poly5/PTLens), TCA (linear/poly3) and vignetting (pa) profiles, lensfun-format loader
Profiles are interpolated in focal length (vignetting also in aperture and focus distance), scaled by crop factor, and converted by least squares to the DNG WarpRectilinear / FixVignetteRadial polynomials the pipeline already applies. The lensfun database (CC-BY-SA) is only read from a directory the user supplies; nothing is bundled. Sources: lensfun's published database format and lens calibration documentation (no lensfun code read), PanoTools wiki (a, b, c distortion model), DNG specification 1.6 (WarpRectilinear, FixVignetteRadial); interpolation and fitting: own design. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
1 parent
0ba118de92
commit
00ea386f7f
8 files changed
+1056
No files matched your search
Generated
+9
@@ -3291,6 +3291,15 @@ dependencies = [
|
||||
"rayon",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "lightcraft-lensdb"
|
||||
version = "0.5.0"
|
||||
dependencies = [
|
||||
"lightcraft-develop",
|
||||
"lightcraft-geom",
|
||||
"quick-xml 0.38.4",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "lightcraft-mcp"
|
||||
version = "0.6.0"
|
||||
|
||||
@@ -24,6 +24,7 @@ lightcraft-raw = { path = "crates/raw" }
|
||||
lightcraft-codecs = { path = "crates/codecs" }
|
||||
lightcraft-meta = { path = "crates/meta" }
|
||||
lightcraft-develop = { path = "crates/develop" }
|
||||
lightcraft-lensdb = { path = "crates/lensdb" }
|
||||
lightcraft-scenes = { path = "crates/scenes" }
|
||||
lightcraft-pipeline = { path = "crates/pipeline" }
|
||||
lightcraft-gpu = { path = "crates/gpu" }
|
||||
|
||||
@@ -0,0 +1,14 @@
|
||||
[package]
|
||||
name = "lightcraft-lensdb"
|
||||
version.workspace = true
|
||||
edition.workspace = true
|
||||
license.workspace = true
|
||||
description = "Lens correction profiles: distortion, TCA and vignetting models, lensfun-format database loader (L2)"
|
||||
|
||||
[lints]
|
||||
workspace = true
|
||||
|
||||
[dependencies]
|
||||
lightcraft-develop = { workspace = true }
|
||||
lightcraft-geom = { workspace = true }
|
||||
quick-xml = { workspace = true }
|
||||
@@ -0,0 +1,458 @@
|
||||
//! Lens correction profiles: distortion, transverse chromatic aberration (TCA) and vignetting.
|
||||
//!
|
||||
//! - [`model`]: the radial models (distortion poly3 / poly5 / PTLens, TCA linear / poly3, vignetting `pa`), and the
|
||||
//! least-squares conversion to the DNG opcode polynomials (`WarpRectilinear`, `FixVignetteRadial`) the pipeline
|
||||
//! already applies in its single geometry resample.
|
||||
//! - [`LensDb`]: profiles read from lensfun-format XML files in a directory the user supplies. The lensfun database
|
||||
//! is CC-BY-SA; it is never bundled with the code or the binary, only read at run time from where the user put it.
|
||||
//! - [`LensDb::find_camera`] / [`LensDb::find_lens`]: identification from EXIF make / model / lens model strings,
|
||||
//! with mount and crop-factor handling.
|
||||
//! - [`Profile::correction`]: the profile at a shot's focal length (and aperture / distance for vignetting),
|
||||
//! interpolated between calibrated points, expressed as [`lightcraft_develop::EmbeddedLens`] for a given image.
|
||||
//!
|
||||
//! Sources: the lensfun database format and model formulas as described in lensfun's published documentation
|
||||
//! (database file format and lens calibration pages; no lensfun code was read), Panorama Tools' distortion model
|
||||
//! (a, b, c) as documented by the PanoTools wiki, Adobe DNG specification 1.6 (opcodes `WarpRectilinear`,
|
||||
//! `FixVignetteRadial`). Interpolation and fitting: own design.
|
||||
#![forbid(unsafe_code)]
|
||||
#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
|
||||
|
||||
pub mod model;
|
||||
mod xml;
|
||||
|
||||
use std::path::Path;
|
||||
|
||||
use lightcraft_develop::{EmbeddedLens, EmbeddedVignette, EmbeddedWarp};
|
||||
use lightcraft_geom::Point;
|
||||
pub use model::{Distortion, Tca, Vignetting};
|
||||
|
||||
/// Largest database file read (lensfun's biggest file is well under 1 MB).
|
||||
const MAX_FILE_BYTES: u64 = 16 << 20;
|
||||
/// Most files read from one directory.
|
||||
const MAX_FILES: usize = 2000;
|
||||
/// Largest radial error (half-diagonal units) accepted when converting a distortion to the DNG polynomial.
|
||||
pub const MAX_WARP_ERROR: f64 = 0.001;
|
||||
/// Largest relative error accepted when converting vignetting to the DNG polynomial.
|
||||
pub const MAX_VIGNETTE_ERROR: f64 = 0.01;
|
||||
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub enum Error {
|
||||
Io(String),
|
||||
Xml(String),
|
||||
TooLarge,
|
||||
}
|
||||
|
||||
impl std::fmt::Display for Error {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
match self {
|
||||
Error::Io(e) => write!(f, "lens database: {e}"),
|
||||
Error::Xml(e) => write!(f, "lens database: invalid XML: {e}"),
|
||||
Error::TooLarge => write!(f, "lens database: file too large or too deeply nested"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl std::error::Error for Error {}
|
||||
|
||||
/// A camera body: what is needed to pick lenses (mount) and scale profiles (crop factor).
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct Camera {
|
||||
pub maker: String,
|
||||
pub model: String,
|
||||
pub mount: String,
|
||||
pub crop_factor: f64,
|
||||
}
|
||||
|
||||
/// One calibrated vignetting point.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct VignettingPoint {
|
||||
pub focal: f64,
|
||||
pub aperture: f64,
|
||||
/// Focus distance, metres (lensfun uses 1000 for infinity).
|
||||
pub distance: f64,
|
||||
pub model: Vignetting,
|
||||
}
|
||||
|
||||
/// A lens profile: names, mounts, the calibration frame and the calibrated points.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct Profile {
|
||||
pub maker: String,
|
||||
/// Every name the lens goes by (EXIF variants, translations).
|
||||
pub models: Vec<String>,
|
||||
pub mounts: Vec<String>,
|
||||
/// Crop factor of the camera the calibration was made on.
|
||||
pub crop_factor: f64,
|
||||
/// Long / short side of the calibration frame.
|
||||
pub aspect_ratio: f64,
|
||||
/// False for fisheye and other non-rectilinear projections: distortion is then not applied.
|
||||
pub rectilinear: bool,
|
||||
/// Sorted by focal length.
|
||||
pub distortion: Vec<(f64, Distortion)>,
|
||||
pub tca: Vec<(f64, Tca)>,
|
||||
pub vignetting: Vec<VignettingPoint>,
|
||||
}
|
||||
|
||||
/// The shot a correction is wanted for.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Shot {
|
||||
/// Focal length, mm.
|
||||
pub focal: f64,
|
||||
pub aperture: Option<f64>,
|
||||
/// Focus distance, metres.
|
||||
pub distance: Option<f64>,
|
||||
/// Crop factor of the camera (from [`LensDb::find_camera`]); the profile's own when unknown.
|
||||
pub crop_factor: Option<f64>,
|
||||
/// Image size in pixels (the frame the correction is expressed for; the optical centre is its centre).
|
||||
pub width: f64,
|
||||
pub height: f64,
|
||||
}
|
||||
|
||||
/// A set of lens and camera profiles.
|
||||
#[derive(Clone, Debug, Default, PartialEq)]
|
||||
pub struct LensDb {
|
||||
pub cameras: Vec<Camera>,
|
||||
pub lenses: Vec<Profile>,
|
||||
}
|
||||
|
||||
impl LensDb {
|
||||
/// Read every `*.xml` file in `dir` (lensfun database format). Unreadable or invalid files are skipped and
|
||||
/// reported in the second value; only a directory that can't be listed is an error.
|
||||
pub fn load_dir(dir: &Path) -> Result<(LensDb, Vec<(String, Error)>), Error> {
|
||||
let mut entries: Vec<_> = std::fs::read_dir(dir)
|
||||
.map_err(|e| Error::Io(format!("{}: {e}", dir.display())))?
|
||||
.flatten()
|
||||
.map(|e| e.path())
|
||||
.filter(|p| p.extension().is_some_and(|x| x.eq_ignore_ascii_case("xml")))
|
||||
.collect();
|
||||
entries.sort();
|
||||
entries.truncate(MAX_FILES);
|
||||
let mut db = LensDb::default();
|
||||
let mut errors = Vec::new();
|
||||
for path in entries {
|
||||
let name = path.display().to_string();
|
||||
match read_capped(&path).and_then(|s| LensDb::parse(&s)) {
|
||||
Ok(part) => db.extend(part),
|
||||
Err(e) => errors.push((name, e)),
|
||||
}
|
||||
}
|
||||
Ok((db, errors))
|
||||
}
|
||||
|
||||
/// Parse one lensfun-format XML document.
|
||||
pub fn parse(xml: &str) -> Result<LensDb, Error> {
|
||||
let root = xml::parse(xml)?;
|
||||
let Some(top) = root.kid("lensdatabase") else { return Err(Error::Xml("no <lensdatabase> element".into())) };
|
||||
let cameras = top.kids("camera").filter_map(camera).collect();
|
||||
let lenses = top.kids("lens").flat_map(lens).collect();
|
||||
Ok(LensDb { cameras, lenses })
|
||||
}
|
||||
|
||||
pub fn extend(&mut self, other: LensDb) {
|
||||
self.cameras.extend(other.cameras);
|
||||
self.lenses.extend(other.lenses);
|
||||
}
|
||||
|
||||
pub fn is_empty(&self) -> bool {
|
||||
self.cameras.is_empty() && self.lenses.is_empty()
|
||||
}
|
||||
|
||||
/// The camera with this EXIF make and model.
|
||||
pub fn find_camera(&self, make: &str, model: &str) -> Option<&Camera> {
|
||||
let (mk, md) = (tokens(make), tokens(model));
|
||||
if md.is_empty() {
|
||||
return None;
|
||||
}
|
||||
self.cameras.iter().find(|c| tokens(&c.model) == md && maker_matches(&tokens(&c.maker), &mk))
|
||||
}
|
||||
|
||||
/// The best profile for an EXIF lens model (and lens make, when the file has one), preferring lenses for the
|
||||
/// camera's mount. `None` when no name is a convincing match.
|
||||
pub fn find_lens(&self, lens_make: Option<&str>, lens_model: &str, camera: Option<&Camera>) -> Option<&Profile> {
|
||||
let exif = tokens(lens_model);
|
||||
if exif.is_empty() {
|
||||
return None;
|
||||
}
|
||||
let make = lens_make.map(tokens).unwrap_or_default();
|
||||
let mut best: Option<(f64, &Profile)> = None;
|
||||
for lens in &self.lenses {
|
||||
let maker = tokens(&lens.maker);
|
||||
if !make.is_empty() && !maker.is_empty() && !maker_matches(&maker, &make) {
|
||||
continue;
|
||||
}
|
||||
let mut score = lens.models.iter().map(|m| name_score(&tokens(m), &exif, &maker)).fold(0.0, f64::max);
|
||||
if let Some(cam) = camera
|
||||
&& !lens.mounts.is_empty()
|
||||
&& !lens.mounts.iter().any(|m| m.eq_ignore_ascii_case(&cam.mount))
|
||||
{
|
||||
score *= 0.8;
|
||||
}
|
||||
if score >= 0.6 && best.is_none_or(|(s, _)| score > s) {
|
||||
best = Some((score, lens));
|
||||
}
|
||||
}
|
||||
best.map(|(_, l)| l)
|
||||
}
|
||||
}
|
||||
|
||||
fn read_capped(path: &Path) -> Result<String, Error> {
|
||||
let len = std::fs::metadata(path).map_err(|e| Error::Io(e.to_string()))?.len();
|
||||
if len > MAX_FILE_BYTES {
|
||||
return Err(Error::TooLarge);
|
||||
}
|
||||
let bytes = std::fs::read(path).map_err(|e| Error::Io(e.to_string()))?;
|
||||
Ok(String::from_utf8_lossy(&bytes).into_owned())
|
||||
}
|
||||
|
||||
/// Lower-case word tokens; letters and digits are split ("24-105mm" → 24, 105, mm), decimals kept ("2.8").
|
||||
fn tokens(s: &str) -> Vec<String> {
|
||||
let mut out: Vec<String> = Vec::new();
|
||||
let mut cur = String::new();
|
||||
let mut prev_digit: Option<bool> = None;
|
||||
for c in s.chars().flat_map(char::to_lowercase) {
|
||||
let digit = c.is_ascii_digit() || (c == '.' && prev_digit == Some(true));
|
||||
if c.is_alphanumeric() || digit {
|
||||
if prev_digit.is_some_and(|p| p != digit) && !cur.is_empty() {
|
||||
out.push(std::mem::take(&mut cur));
|
||||
}
|
||||
cur.push(c);
|
||||
prev_digit = Some(digit);
|
||||
} else {
|
||||
if !cur.is_empty() {
|
||||
out.push(std::mem::take(&mut cur));
|
||||
}
|
||||
prev_digit = None;
|
||||
}
|
||||
}
|
||||
if !cur.is_empty() {
|
||||
out.push(cur);
|
||||
}
|
||||
for t in &mut out {
|
||||
// "4.0" and "4" are the same aperture
|
||||
if t.contains('.') && t.starts_with(|c: char| c.is_ascii_digit()) {
|
||||
let trimmed = t.trim_end_matches('0').trim_end_matches('.');
|
||||
*t = trimmed.to_string();
|
||||
}
|
||||
}
|
||||
out.retain(|t| !t.is_empty());
|
||||
out
|
||||
}
|
||||
|
||||
fn is_number(t: &str) -> bool {
|
||||
t.starts_with(|c: char| c.is_ascii_digit())
|
||||
}
|
||||
|
||||
/// Makers match when one's first word appears in the other ("NIKON CORPORATION" / "Nikon").
|
||||
fn maker_matches(a: &[String], b: &[String]) -> bool {
|
||||
match (a.first(), b.first()) {
|
||||
(Some(x), Some(y)) => b.contains(x) || a.contains(y),
|
||||
_ => true,
|
||||
}
|
||||
}
|
||||
|
||||
/// How well a database name matches the EXIF name, 0..1. Every number in the database name (focal lengths,
|
||||
/// apertures) must appear in the EXIF name; the maker's words are ignored on both sides (EXIF often omits them).
|
||||
fn name_score(db: &[String], exif: &[String], maker: &[String]) -> f64 {
|
||||
let db: Vec<&String> = db.iter().filter(|t| !maker.contains(t)).collect();
|
||||
let exif: Vec<&String> = exif.iter().filter(|t| !maker.contains(t)).collect();
|
||||
if db.is_empty() || exif.is_empty() {
|
||||
return 0.0;
|
||||
}
|
||||
if db.iter().any(|t| is_number(t) && !exif.contains(t)) || exif.iter().any(|t| is_number(t) && !db.contains(t)) {
|
||||
return 0.0;
|
||||
}
|
||||
let common = db.iter().filter(|t| exif.contains(t)).count() as f64;
|
||||
common / db.len().max(exif.len()) as f64
|
||||
}
|
||||
|
||||
fn camera(n: &xml::Node) -> Option<Camera> {
|
||||
Some(Camera {
|
||||
maker: n.kid_text("maker")?.to_string(),
|
||||
model: n.kid_text("model")?.to_string(),
|
||||
mount: n.kid_text("mount").unwrap_or_default().to_string(),
|
||||
crop_factor: crop_of(n.kid_text("cropfactor")).unwrap_or(1.0),
|
||||
})
|
||||
}
|
||||
|
||||
fn crop_of(s: Option<&str>) -> Option<f64> {
|
||||
s.and_then(|v| v.trim().parse::<f64>().ok()).filter(|v| v.is_finite() && (0.05..=20.0).contains(v))
|
||||
}
|
||||
|
||||
fn aspect_of(s: Option<&str>) -> Option<f64> {
|
||||
let s = s?.trim();
|
||||
let v = match s.split_once(':') {
|
||||
Some((a, b)) => a.trim().parse::<f64>().ok()? / b.trim().parse::<f64>().ok()?,
|
||||
None => s.parse::<f64>().ok()?,
|
||||
};
|
||||
let v = if v < 1.0 { 1.0 / v } else { v };
|
||||
(v.is_finite() && (1.0..=10.0).contains(&v)).then_some(v)
|
||||
}
|
||||
|
||||
fn positive(v: Option<f64>) -> Option<f64> {
|
||||
v.filter(|v| *v > 0.0 && *v < 1e6)
|
||||
}
|
||||
|
||||
/// One `<lens>` element → one profile per `<calibration>` (newer files give each its own crop factor and aspect).
|
||||
fn lens(n: &xml::Node) -> Vec<Profile> {
|
||||
let Some(maker) = n.kid_text("maker") else { return Vec::new() };
|
||||
let models: Vec<String> = n.kids("model").map(|m| m.text.trim().to_string()).filter(|m| !m.is_empty()).collect();
|
||||
if models.is_empty() {
|
||||
return Vec::new();
|
||||
}
|
||||
let crop = crop_of(n.kid_text("cropfactor")).unwrap_or(1.0);
|
||||
let aspect = aspect_of(n.kid_text("aspect-ratio")).unwrap_or(1.5);
|
||||
let rectilinear = n.kid_text("type").is_none_or(|t| t.eq_ignore_ascii_case("rectilinear"));
|
||||
let mounts = n.kids("mount").map(|m| m.text.trim().to_string()).filter(|m| !m.is_empty()).collect::<Vec<_>>();
|
||||
n.kids("calibration")
|
||||
.map(|cal| {
|
||||
let mut p = Profile {
|
||||
maker: maker.to_string(),
|
||||
models: models.clone(),
|
||||
mounts: mounts.clone(),
|
||||
crop_factor: crop_of(cal.attr("cropfactor")).unwrap_or(crop),
|
||||
aspect_ratio: aspect_of(cal.attr("aspect-ratio")).unwrap_or(aspect),
|
||||
rectilinear,
|
||||
distortion: Vec::new(),
|
||||
tca: Vec::new(),
|
||||
vignetting: Vec::new(),
|
||||
};
|
||||
for c in &cal.children {
|
||||
let Some(focal) = positive(c.num("focal")) else { continue };
|
||||
let k = |name: &str| c.num(name).unwrap_or(0.0);
|
||||
match (c.name.as_str(), c.attr("model")) {
|
||||
("distortion", Some("poly3")) => p.distortion.push((focal, Distortion::Poly3 { k1: k("k1") })),
|
||||
("distortion", Some("poly5")) => p.distortion.push((focal, Distortion::Poly5 { k1: k("k1"), k2: k("k2") })),
|
||||
("distortion", Some("ptlens")) => p.distortion.push((focal, Distortion::PtLens { a: k("a"), b: k("b"), c: k("c") })),
|
||||
("tca", Some("linear")) => {
|
||||
let one = |name: &str| c.num(name).unwrap_or(1.0);
|
||||
p.tca.push((focal, Tca::Linear { kr: one("kr"), kb: one("kb") }));
|
||||
}
|
||||
("tca", Some("poly3")) => {
|
||||
let one = |name: &str| c.num(name).unwrap_or(1.0);
|
||||
p.tca.push((focal, Tca::Poly3 { vr: one("vr"), vb: one("vb"), cr: k("cr"), cb: k("cb"), br: k("br"), bb: k("bb") }));
|
||||
}
|
||||
("vignetting", Some("pa")) => {
|
||||
let (Some(aperture), distance) = (positive(c.num("aperture")), positive(c.num("distance")).unwrap_or(1000.0)) else {
|
||||
continue;
|
||||
};
|
||||
let model = Vignetting { k1: k("k1"), k2: k("k2"), k3: k("k3") };
|
||||
p.vignetting.push(VignettingPoint { focal, aperture, distance, model });
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
p.distortion.sort_by(|a, b| a.0.total_cmp(&b.0));
|
||||
p.tca.sort_by(|a, b| a.0.total_cmp(&b.0));
|
||||
p
|
||||
})
|
||||
.filter(|p| !p.distortion.is_empty() || !p.tca.is_empty() || !p.vignetting.is_empty())
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Linear interpolation weights in a sorted key list: the bracketing indices and the weight of the upper one.
|
||||
fn bracket(keys: &[f64], x: f64) -> Option<(usize, usize, f64)> {
|
||||
let last = keys.len().checked_sub(1)?;
|
||||
let hi = keys.iter().position(|&k| k >= x).unwrap_or(last);
|
||||
let lo = hi.saturating_sub(1);
|
||||
let (a, b) = (keys.get(lo).copied()?, keys.get(hi).copied()?);
|
||||
if x <= a || b <= a {
|
||||
return Some((lo, lo, 0.0));
|
||||
}
|
||||
if x >= b {
|
||||
return Some((hi, hi, 0.0));
|
||||
}
|
||||
Some((lo, hi, (x - a) / (b - a)))
|
||||
}
|
||||
|
||||
/// A value interpolated in focal length between calibrated points, by evaluating both neighbours.
|
||||
fn by_focal<T: Copy>(points: &[(f64, T)], focal: f64, eval: impl Fn(&T) -> f64) -> Option<f64> {
|
||||
let keys: Vec<f64> = points.iter().map(|p| p.0).collect();
|
||||
let (lo, hi, t) = bracket(&keys, focal)?;
|
||||
let (a, b) = (eval(&points.get(lo)?.1), eval(&points.get(hi)?.1));
|
||||
Some(a + (b - a) * t)
|
||||
}
|
||||
|
||||
impl Profile {
|
||||
/// Distortion ratio `r_d / r_u` at lensfun radius `r` for `focal` (1 when the profile has none).
|
||||
pub fn distortion_ratio(&self, focal: f64, r: f64) -> f64 {
|
||||
if !self.rectilinear {
|
||||
return 1.0;
|
||||
}
|
||||
by_focal(&self.distortion, focal, |d| d.ratio(r)).unwrap_or(1.0)
|
||||
}
|
||||
|
||||
/// TCA ratio for channel `ch` at lensfun radius `r` for `focal`.
|
||||
pub fn tca_ratio(&self, focal: f64, ch: usize, r: f64) -> f64 {
|
||||
by_focal(&self.tca, focal, |t| t.ratio(ch, r)).unwrap_or(1.0)
|
||||
}
|
||||
|
||||
/// Vignetting attenuation at radius `r` (half-diagonal units): linear in focal length, in stops of aperture and
|
||||
/// in 1 / distance between the calibrated points.
|
||||
pub fn vignetting_attenuation(&self, focal: f64, aperture: Option<f64>, distance: Option<f64>, r: f64) -> Option<f64> {
|
||||
let mut focals: Vec<f64> = self.vignetting.iter().map(|v| v.focal).collect();
|
||||
focals.sort_by(f64::total_cmp);
|
||||
focals.dedup();
|
||||
let (lo, hi, t) = bracket(&focals, focal)?;
|
||||
let at = |f: f64| -> Option<f64> {
|
||||
let pts: Vec<&VignettingPoint> = self.vignetting.iter().filter(|v| v.focal == f).collect();
|
||||
// the calibrated distance nearest the shot's (in 1/d; infinity when unknown)
|
||||
let inv = |d: f64| 1.0 / d.max(0.01);
|
||||
let want = distance.map_or(0.0, inv);
|
||||
let near = pts.iter().map(|v| v.distance).min_by(|a, b| (inv(*a) - want).abs().total_cmp(&(inv(*b) - want).abs()))?;
|
||||
let mut at_d: Vec<&&VignettingPoint> = pts.iter().filter(|v| v.distance == near).collect();
|
||||
at_d.sort_by(|a, b| a.aperture.total_cmp(&b.aperture));
|
||||
let stops: Vec<f64> = at_d.iter().map(|v| v.aperture.log2()).collect();
|
||||
// unknown aperture: the widest calibrated one (the strongest falloff)
|
||||
let n = aperture.filter(|a| *a > 0.0).map_or(stops.first().copied()?, f64::log2);
|
||||
let (i, j, u) = bracket(&stops, n)?;
|
||||
let (a, b) = (at_d.get(i)?.model.attenuation(r), at_d.get(j)?.model.attenuation(r));
|
||||
Some(a + (b - a) * u)
|
||||
};
|
||||
let (a, b) = (at(*focals.get(lo)?)?, at(*focals.get(hi)?)?);
|
||||
Some(a + (b - a) * t)
|
||||
}
|
||||
|
||||
/// The correction for `shot`, as the DNG opcode polynomials the pipeline applies, centred on the image.
|
||||
/// Distortion and TCA are combined per plane (R, G, B); a model the polynomial can't follow closely
|
||||
/// ([`MAX_WARP_ERROR`], [`MAX_VIGNETTE_ERROR`]) or that would fold the image is left out.
|
||||
pub fn correction(&self, shot: &Shot) -> Option<EmbeddedLens> {
|
||||
let (w, h) = (shot.width, shot.height);
|
||||
if !(w >= 2.0 && h >= 2.0 && w.is_finite() && h.is_finite() && shot.focal.is_finite() && shot.focal > 0.0) {
|
||||
return None;
|
||||
}
|
||||
let cam = shot.crop_factor.filter(|c| c.is_finite() && *c > 0.0).unwrap_or(self.crop_factor);
|
||||
// DNG radius (1 = image half diagonal) → lensfun radii
|
||||
let to_vig = self.crop_factor / cam;
|
||||
let to_dist = to_vig * (1.0 + self.aspect_ratio * self.aspect_ratio).sqrt();
|
||||
let center = Point::new(0.5, 0.5);
|
||||
let radius = w.hypot(h) / 2.0 / w.max(h);
|
||||
let f = shot.focal;
|
||||
let mut lens = EmbeddedLens::default();
|
||||
let has_dist = self.rectilinear && !self.distortion.is_empty();
|
||||
if has_dist || !self.tca.is_empty() {
|
||||
let plane = |ch: usize| -> Option<[f64; 6]> {
|
||||
let ratio = |r: f64| {
|
||||
let rl = r * to_dist;
|
||||
let d = self.distortion_ratio(f, rl);
|
||||
d * self.tca_ratio(f, ch, rl * d)
|
||||
};
|
||||
let (k, err) = model::fit_even4(ratio)?;
|
||||
(err <= MAX_WARP_ERROR).then_some([k[0], k[1], k[2], k[3], 0.0, 0.0])
|
||||
};
|
||||
if let (Some(r), Some(g), Some(b)) = (plane(0), plane(1), plane(2)) {
|
||||
lens.warp = Some(EmbeddedWarp { planes: [r, g, b], center, radius });
|
||||
}
|
||||
}
|
||||
if !self.vignetting.is_empty() {
|
||||
let att = |r: f64| self.vignetting_attenuation(f, shot.aperture, shot.distance, r * to_vig).unwrap_or(1.0);
|
||||
if let Some((k, err)) = model::fit_gain5(|r| 1.0 / att(r).max(0.05))
|
||||
&& err <= MAX_VIGNETTE_ERROR
|
||||
{
|
||||
lens.vignette = Some(EmbeddedVignette { k, center, radius });
|
||||
}
|
||||
}
|
||||
(lens.warp.is_some() || lens.vignette.is_some()).then_some(lens)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
@@ -0,0 +1,192 @@
|
||||
//! The correction models and their conversion to the DNG-style polynomials the pipeline applies.
|
||||
//!
|
||||
//! All models are functions of a normalized radius `r`. Following the lensfun database documentation, distortion and
|
||||
//! TCA use `r = 1` at half the shorter side of the calibration frame, vignetting `r = 1` at half its diagonal (the
|
||||
//! corner). Distortion and TCA map a corrected (undistorted) radius to the radius where the light actually landed,
|
||||
//! which is the inverse map a resampler needs.
|
||||
|
||||
/// Radial distortion: `r_d = r_u · ratio(r_u)`.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub enum Distortion {
|
||||
/// `r_d = r_u · (1 − k1 + k1·r_u²)`.
|
||||
Poly3 { k1: f64 },
|
||||
/// `r_d = r_u · (1 + k1·r_u² + k2·r_u⁴)`.
|
||||
Poly5 { k1: f64, k2: f64 },
|
||||
/// Panorama Tools: `r_d = r_u · (a·r_u³ + b·r_u² + c·r_u + 1 − a − b − c)`.
|
||||
PtLens { a: f64, b: f64, c: f64 },
|
||||
}
|
||||
|
||||
impl Distortion {
|
||||
/// `r_d / r_u` at `r_u`.
|
||||
pub fn ratio(&self, r: f64) -> f64 {
|
||||
match *self {
|
||||
Distortion::Poly3 { k1 } => 1.0 - k1 + k1 * r * r,
|
||||
Distortion::Poly5 { k1, k2 } => {
|
||||
let r2 = r * r;
|
||||
1.0 + k1 * r2 + k2 * r2 * r2
|
||||
}
|
||||
Distortion::PtLens { a, b, c } => ((a * r + b) * r + c) * r + 1.0 - a - b - c,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Transverse (lateral) chromatic aberration of red and blue relative to green: `r_ch = r_g · ratio(ch, r_g)`.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub enum Tca {
|
||||
/// `r_ch = r_g · k`.
|
||||
Linear { kr: f64, kb: f64 },
|
||||
/// `r_ch = r_g · (b·r_g² + c·r_g + v)`, per channel.
|
||||
Poly3 { vr: f64, vb: f64, cr: f64, cb: f64, br: f64, bb: f64 },
|
||||
}
|
||||
|
||||
impl Tca {
|
||||
/// `r_ch / r_g` for `ch` 0 (red), 2 (blue); 1 for green.
|
||||
pub fn ratio(&self, ch: usize, r: f64) -> f64 {
|
||||
match (*self, ch) {
|
||||
(_, 1) => 1.0,
|
||||
(Tca::Linear { kr, .. }, 0) => kr,
|
||||
(Tca::Linear { kb, .. }, _) => kb,
|
||||
(Tca::Poly3 { vr, cr, br, .. }, 0) => (br * r + cr) * r + vr,
|
||||
(Tca::Poly3 { vb, cb, bb, .. }, _) => (bb * r + cb) * r + vb,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Vignetting, `pa` model: the recorded brightness is the true brightness times `1 + k1·r² + k2·r⁴ + k3·r⁶`.
|
||||
#[derive(Clone, Copy, Debug, Default, PartialEq)]
|
||||
pub struct Vignetting {
|
||||
pub k1: f64,
|
||||
pub k2: f64,
|
||||
pub k3: f64,
|
||||
}
|
||||
|
||||
impl Vignetting {
|
||||
pub fn attenuation(&self, r: f64) -> f64 {
|
||||
let r2 = r * r;
|
||||
1.0 + r2 * (self.k1 + r2 * (self.k2 + r2 * self.k3))
|
||||
}
|
||||
}
|
||||
|
||||
/// Least squares `Σ w (Σ_j c_j·basis_j(x) − y)²` for small dense systems (normal equations, partial pivoting).
|
||||
/// `None` when the system is singular or a value isn't finite.
|
||||
pub(crate) fn least_squares<const N: usize>(samples: &[([f64; N], f64)]) -> Option<[f64; N]> {
|
||||
let mut a = [[0.0; N]; N];
|
||||
let mut b = [0.0; N];
|
||||
for (basis, y) in samples {
|
||||
for i in 0..N {
|
||||
b[i] += basis[i] * y;
|
||||
for j in 0..N {
|
||||
a[i][j] += basis[i] * basis[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
for col in 0..N {
|
||||
let pivot = (col..N).max_by(|&x, &y| a[x][col].abs().total_cmp(&a[y][col].abs()))?;
|
||||
if !a[pivot][col].is_finite() || a[pivot][col].abs() < 1e-14 {
|
||||
return None;
|
||||
}
|
||||
a.swap(col, pivot);
|
||||
b.swap(col, pivot);
|
||||
for row in 0..N {
|
||||
if row != col {
|
||||
let f = a[row][col] / a[col][col];
|
||||
for k in col..N {
|
||||
a[row][k] -= f * a[col][k];
|
||||
}
|
||||
b[row] -= f * b[col];
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut x = [0.0; N];
|
||||
for i in 0..N {
|
||||
x[i] = b[i] / a[i][i];
|
||||
}
|
||||
x.iter().all(|v| v.is_finite()).then_some(x)
|
||||
}
|
||||
|
||||
/// Fit `ratio(r) ≈ k0 + k1·r² + k2·r⁴ + k3·r⁶` over `r ∈ [0, 1]` (the DNG `WarpRectilinear` radial form).
|
||||
/// Returns the coefficients and the largest radial error `max |r·(fit − ratio)|`, in the same units as `r`.
|
||||
pub fn fit_even4(ratio: impl Fn(f64) -> f64) -> Option<([f64; 4], f64)> {
|
||||
let pts: Vec<f64> = (0..=SAMPLES).map(|i| i as f64 / SAMPLES as f64).collect();
|
||||
let samples: Vec<([f64; 4], f64)> = pts
|
||||
.iter()
|
||||
.map(|&r| {
|
||||
let r2 = r * r;
|
||||
// weight by r: the error that matters is the displacement r·Δratio
|
||||
([r, r * r2, r * r2 * r2, r * r2 * r2 * r2], r * ratio(r))
|
||||
})
|
||||
.collect();
|
||||
let k = least_squares(&samples)?;
|
||||
let eval = |r: f64| {
|
||||
let r2 = r * r;
|
||||
k[0] + r2 * (k[1] + r2 * (k[2] + r2 * k[3]))
|
||||
};
|
||||
let err = pts.iter().map(|&r| (r * (eval(r) - ratio(r))).abs()).fold(0.0, f64::max);
|
||||
// the warp must stay monotonic (no folding): r·fit(r) increasing
|
||||
let mono = pts.windows(2).all(|w| w[1] * eval(w[1]) > w[0] * eval(w[0]));
|
||||
(err.is_finite() && mono).then_some((k, err))
|
||||
}
|
||||
|
||||
/// Fit `gain(r) ≈ 1 + k0·r² + … + k4·r¹⁰` over `r ∈ [0, 1]` (the DNG `FixVignetteRadial` form).
|
||||
/// Returns the coefficients and the largest relative gain error.
|
||||
pub fn fit_gain5(gain: impl Fn(f64) -> f64) -> Option<([f64; 5], f64)> {
|
||||
let pts: Vec<f64> = (0..=SAMPLES).map(|i| i as f64 / SAMPLES as f64).collect();
|
||||
let samples: Vec<([f64; 5], f64)> = pts
|
||||
.iter()
|
||||
.map(|&r| {
|
||||
let r2 = r * r;
|
||||
([r2, r2 * r2, r2.powi(3), r2.powi(4), r2.powi(5)], gain(r) - 1.0)
|
||||
})
|
||||
.collect();
|
||||
let k = least_squares(&samples)?;
|
||||
let eval = |r: f64| {
|
||||
let r2 = r * r;
|
||||
1.0 + r2 * (k[0] + r2 * (k[1] + r2 * (k[2] + r2 * (k[3] + r2 * k[4]))))
|
||||
};
|
||||
let err = pts.iter().map(|&r| ((eval(r) - gain(r)) / gain(r).abs().max(1e-9)).abs()).fold(0.0, f64::max);
|
||||
err.is_finite().then_some((k, err))
|
||||
}
|
||||
|
||||
const SAMPLES: usize = 64;
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn models_are_identity_at_centre_and_documented_at_one() {
|
||||
assert_eq!(Distortion::Poly3 { k1: 0.1 }.ratio(1.0), 1.0, "poly3 keeps r = 1");
|
||||
assert!((Distortion::Poly5 { k1: 0.1, k2: 0.01 }.ratio(1.0) - 1.11).abs() < 1e-12);
|
||||
assert!((Distortion::PtLens { a: 0.01, b: -0.02, c: 0.03 }.ratio(1.0) - 1.0).abs() < 1e-12, "ptlens keeps r = 1");
|
||||
assert_eq!(Tca::Linear { kr: 1.001, kb: 0.999 }.ratio(0, 0.5), 1.001);
|
||||
assert_eq!(Tca::Poly3 { vr: 1.0, vb: 1.0, cr: 0.0, cb: 0.0, br: 0.001, bb: 0.0 }.ratio(1, 0.7), 1.0);
|
||||
assert!((Vignetting { k1: -0.3, k2: 0.0, k3: 0.0 }.attenuation(1.0) - 0.7).abs() < 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn even_fit_is_exact_for_even_models_and_close_for_ptlens() {
|
||||
let (k, err) = fit_even4(|r| Distortion::Poly5 { k1: -0.05, k2: 0.01 }.ratio(r)).unwrap();
|
||||
assert!(err < 1e-9, "{err}");
|
||||
assert!((k[0] - 1.0).abs() < 1e-9 && (k[1] + 0.05).abs() < 1e-9 && (k[2] - 0.01).abs() < 1e-9);
|
||||
let (_, err) = fit_even4(|r| Distortion::PtLens { a: 0.01, b: -0.03, c: 0.0 }.ratio(r)).unwrap();
|
||||
assert!(err < 2e-4, "{err}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn folding_maps_are_refused() {
|
||||
assert!(fit_even4(|r| 1.0 - 2.0 * r * r).is_none());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn gain_fit_matches_inverse_attenuation() {
|
||||
let v = Vignetting { k1: -0.4, k2: 0.1, k3: -0.02 };
|
||||
let (_, err) = fit_gain5(|r| 1.0 / v.attenuation(r)).unwrap();
|
||||
assert!(err < 1e-3, "{err}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn singular_systems_are_none() {
|
||||
assert!(least_squares::<2>(&[([1.0, 1.0], 1.0), ([2.0, 2.0], 2.0)]).is_none());
|
||||
assert!(least_squares::<1>(&[]).is_none());
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,245 @@
|
||||
use super::*;
|
||||
|
||||
/// A synthetic database in the lensfun format (made-up coefficients, not lensfun data).
|
||||
const DB: &str = r#"<?xml version="1.0" encoding="UTF-8"?>
|
||||
<lensdatabase version="2">
|
||||
<mount><name>Test Mount</name></mount>
|
||||
<camera>
|
||||
<maker>Testcam Corporation</maker>
|
||||
<model>TC-1</model>
|
||||
<mount>Test Mount</mount>
|
||||
<cropfactor>1.5</cropfactor>
|
||||
</camera>
|
||||
<camera>
|
||||
<maker>Testcam</maker>
|
||||
<model>TC-FF</model>
|
||||
<mount>Test Mount</mount>
|
||||
<cropfactor>1</cropfactor>
|
||||
</camera>
|
||||
<lens>
|
||||
<maker>Testcam</maker>
|
||||
<model>Testcam 24-70mm f/2.8 & Co</model>
|
||||
<model lang="de">Testcam 24-70mm f/2.8 Zoom</model>
|
||||
<mount>Test Mount</mount>
|
||||
<cropfactor>1</cropfactor>
|
||||
<calibration>
|
||||
<distortion model="ptlens" focal="24" a="0.01" b="-0.04" c="0" />
|
||||
<distortion model="poly3" focal="70" k1="0.01" />
|
||||
<tca model="linear" focal="24" kr="1.0004" kb="0.9997" />
|
||||
<tca model="poly3" focal="70" vr="1.0001" vb="0.9999" cr="0" cb="0" br="0.0001" bb="-0.0001" />
|
||||
<vignetting model="pa" focal="24" aperture="2.8" distance="1000" k1="-0.6" k2="0.2" k3="-0.05" />
|
||||
<vignetting model="pa" focal="24" aperture="8" distance="1000" k1="-0.2" k2="0" k3="0" />
|
||||
<vignetting model="pa" focal="70" aperture="2.8" distance="1000" k1="-0.3" k2="0" k3="0" />
|
||||
</calibration>
|
||||
</lens>
|
||||
<lens>
|
||||
<maker>Testcam</maker>
|
||||
<model>Testcam 50mm f/1.8</model>
|
||||
<mount>Other Mount</mount>
|
||||
<calibration>
|
||||
<distortion model="poly5" focal="50" k1="-0.02" k2="0.0003" />
|
||||
</calibration>
|
||||
</lens>
|
||||
<lens>
|
||||
<maker>Fishco</maker>
|
||||
<model>Fishco 8mm f/3.5 Fisheye</model>
|
||||
<type>equisolid</type>
|
||||
<calibration>
|
||||
<distortion model="poly3" focal="8" k1="0.5" />
|
||||
</calibration>
|
||||
</lens>
|
||||
<lens><maker>Broken</maker><model>Broken 1mm</model><calibration><distortion model="poly3" focal="nan" k1="1"/></calibration></lens>
|
||||
</lensdatabase>"#;
|
||||
|
||||
fn db() -> LensDb {
|
||||
LensDb::parse(DB).unwrap()
|
||||
}
|
||||
|
||||
fn shot(focal: f64) -> Shot {
|
||||
Shot { focal, aperture: Some(2.8), distance: None, crop_factor: Some(1.0), width: 6000.0, height: 4000.0 }
|
||||
}
|
||||
|
||||
/// The DNG `WarpRectilinear` source position (the pipeline's formula, radial part).
|
||||
fn dng_source(wp: &EmbeddedWarp, ch: usize, x: f64, y: f64, w: f64, h: f64) -> (f64, f64) {
|
||||
let (cx, cy) = (wp.center.x * w, wp.center.y * h);
|
||||
let m = wp.radius * w.max(h);
|
||||
let (dx, dy) = ((x - cx) / m, (y - cy) / m);
|
||||
let r2 = dx * dx + dy * dy;
|
||||
let [k0, k1, k2, k3, ..] = wp.planes[ch];
|
||||
let f = k0 + r2 * (k1 + r2 * (k2 + r2 * k3));
|
||||
(cx + dx * f * m, cy + dy * f * m)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn parses_cameras_lenses_and_skips_invalid_entries() {
|
||||
let db = db();
|
||||
assert_eq!(db.cameras.len(), 2);
|
||||
assert_eq!(db.lenses.len(), 3, "the lens without valid calibration points is dropped");
|
||||
let z = &db.lenses[0];
|
||||
assert_eq!(z.models[0], "Testcam 24-70mm f/2.8 & Co", "entities decoded");
|
||||
assert_eq!(z.distortion.len(), 2);
|
||||
assert_eq!(z.tca.len(), 2);
|
||||
assert_eq!(z.vignetting.len(), 3);
|
||||
assert_eq!(z.aspect_ratio, 1.5);
|
||||
assert!(!db.lenses[2].rectilinear);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn malformed_documents_are_errors_not_panics() {
|
||||
assert!(LensDb::parse("").is_err());
|
||||
assert!(LensDb::parse("<lensdatabase><lens><maker>").is_ok_and(|d| d.lenses.is_empty()));
|
||||
assert!(LensDb::parse("<other/>").is_err());
|
||||
let deep = "<a>".repeat(100);
|
||||
assert_eq!(LensDb::parse(&deep), Err(Error::TooLarge));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn finds_cameras_by_exif_strings() {
|
||||
let db = db();
|
||||
assert_eq!(db.find_camera("TESTCAM CORPORATION", "TC-1").unwrap().crop_factor, 1.5);
|
||||
assert_eq!(db.find_camera("Testcam", "tc-ff").unwrap().crop_factor, 1.0);
|
||||
assert!(db.find_camera("Testcam", "TC-2").is_none());
|
||||
assert!(db.find_camera("Other", "TC-1").is_none());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn finds_lenses_by_exif_model_with_or_without_the_maker() {
|
||||
let db = db();
|
||||
let cam = db.find_camera("Testcam", "TC-FF");
|
||||
let l = db.find_lens(None, "24-70mm F2.8 Zoom", cam).unwrap();
|
||||
assert_eq!(l.maker, "Testcam");
|
||||
assert!(db.find_lens(Some("Testcam"), "Testcam 50mm F1.8", cam).is_some());
|
||||
assert!(db.find_lens(None, "24-105mm F4", cam).is_none(), "focal lengths must agree");
|
||||
assert!(db.find_lens(Some("Otherco"), "24-70mm f/2.8 Zoom", cam).is_none(), "lens maker must agree");
|
||||
assert!(db.find_lens(None, "", cam).is_none());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn distortion_matches_the_model_within_half_a_pixel_on_a_grid() {
|
||||
let db = db();
|
||||
let lens = &db.lenses[0];
|
||||
let (w, h) = (6000.0, 4000.0);
|
||||
for focal in [24.0, 40.0, 70.0] {
|
||||
let wp = lens.correction(&shot(focal)).unwrap().warp.unwrap();
|
||||
// lensfun normalisation: r = 1 at half the shorter side (full frame, 3:2)
|
||||
let n = h / 2.0;
|
||||
let mut worst: f64 = 0.0;
|
||||
for gy in 0..=20 {
|
||||
for gx in 0..=30 {
|
||||
let (x, y) = (gx as f64 * w / 30.0, gy as f64 * h / 20.0);
|
||||
let (sx, sy) = dng_source(&wp, 1, x, y, w, h);
|
||||
let (dx, dy) = ((x - w / 2.0) / n, (y - h / 2.0) / n);
|
||||
let r = dx.hypot(dy);
|
||||
let k = lens.distortion_ratio(focal, r);
|
||||
let (ex, ey) = (w / 2.0 + dx * k * n, h / 2.0 + dy * k * n);
|
||||
worst = worst.max((sx - ex).hypot(sy - ey));
|
||||
}
|
||||
}
|
||||
assert!(worst < 0.5, "{focal} mm: {worst} px");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn straight_lines_stay_straight_after_correction() {
|
||||
// A synthetic barrel-distorted chart: the camera recorded straight scene lines at r_d = r_u·ratio(r_u).
|
||||
// Correcting means: for each output point of a straight line, the warp fetches the recorded point, which must
|
||||
// be that line's recorded (curved) image. So sample the curved recorded line by inverting the model, then
|
||||
// check the warp sends straight output points onto it.
|
||||
let db = db();
|
||||
let lens = &db.lenses[0];
|
||||
let (w, h) = (6000.0, 4000.0);
|
||||
let wp = lens.correction(&shot(24.0)).unwrap().warp.unwrap();
|
||||
let n = h / 2.0;
|
||||
let y_line = 300.0; // a horizontal line near the top edge
|
||||
let mut dev: f64 = 0.0;
|
||||
for i in 0..=60 {
|
||||
let x = i as f64 * w / 60.0;
|
||||
let (sx, sy) = dng_source(&wp, 1, x, y_line, w, h);
|
||||
// undo the recorded distortion of (sx, sy): find r_u with r_u·ratio(r_u) = r_d (bisection)
|
||||
let (dx, dy) = ((sx - w / 2.0) / n, (sy - h / 2.0) / n);
|
||||
let rd = dx.hypot(dy);
|
||||
let (mut lo, mut hi) = (0.0, 3.0);
|
||||
for _ in 0..80 {
|
||||
let mid = (lo + hi) / 2.0;
|
||||
if mid * lens.distortion_ratio(24.0, mid) < rd {
|
||||
lo = mid;
|
||||
} else {
|
||||
hi = mid;
|
||||
}
|
||||
}
|
||||
let s = if rd > 0.0 { lo / rd } else { 1.0 };
|
||||
let uy = h / 2.0 + dy * s * n;
|
||||
dev = dev.max((uy - y_line).abs());
|
||||
}
|
||||
assert!(dev < 0.5, "line bent by {dev} px");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tca_scales_red_and_blue_planes() {
|
||||
let db = db();
|
||||
let wp = db.lenses[0].correction(&shot(24.0)).unwrap().warp.unwrap();
|
||||
let (w, h) = (6000.0, 4000.0);
|
||||
let (rx, _) = dng_source(&wp, 0, w, h / 2.0, w, h);
|
||||
let (gx, _) = dng_source(&wp, 1, w, h / 2.0, w, h);
|
||||
let (bx, _) = dng_source(&wp, 2, w, h / 2.0, w, h);
|
||||
assert!(((rx - w / 2.0) / (gx - w / 2.0) - 1.0004).abs() < 1e-5, "{rx} {gx}");
|
||||
assert!(((bx - w / 2.0) / (gx - w / 2.0) - 0.9997).abs() < 1e-5, "{bx} {gx}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn vignetting_interpolates_in_aperture_and_focal() {
|
||||
let db = db();
|
||||
let l = &db.lenses[0];
|
||||
let at = |f: f64, n: f64| l.vignetting_attenuation(f, Some(n), None, 1.0).unwrap();
|
||||
assert!((at(24.0, 2.8) - 0.55).abs() < 1e-12);
|
||||
assert!((at(24.0, 8.0) - 0.8).abs() < 1e-12);
|
||||
let mid = at(24.0, 2.8 * 2f64.powf((8f64 / 2.8).log2() / 2.0));
|
||||
assert!((mid - 0.675).abs() < 1e-9, "{mid}");
|
||||
assert!((at(47.0, 2.8) - 0.625).abs() < 1e-9);
|
||||
// the gain undoes the falloff at the corner
|
||||
let v = l.correction(&shot(24.0)).unwrap().vignette.unwrap();
|
||||
let g: f64 = 1.0 + v.k.iter().enumerate().map(|(i, k)| k * 1f64.powi(2 * (i as i32 + 1))).sum::<f64>();
|
||||
assert!((g * 0.55 - 1.0).abs() < 0.01, "{g}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn crop_factor_scales_the_profile() {
|
||||
// On a 1.5× body the image covers the centre of the full-frame calibration: much less distortion at its corner.
|
||||
let db = db();
|
||||
let l = &db.lenses[1];
|
||||
let ff = l.correction(&shot(50.0)).unwrap().warp.unwrap();
|
||||
let aps = l.correction(&Shot { crop_factor: Some(1.5), ..shot(50.0) }).unwrap().warp.unwrap();
|
||||
let corner = |wp: &EmbeddedWarp| dng_source(wp, 1, 0.0, 0.0, 6000.0, 4000.0);
|
||||
let (a, b) = (corner(&ff), corner(&aps));
|
||||
let r = |(x, y): (f64, f64)| (x - 3000.0).hypot(y - 2000.0);
|
||||
let rc = 3000f64.hypot(2000.0);
|
||||
assert!((r(a) - rc).abs() > (r(b) - rc).abs() * 1.5, "{} {}", r(a) - rc, r(b) - rc);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn fisheyes_get_no_rectilinear_warp_and_bad_shots_none() {
|
||||
let db = db();
|
||||
assert!(db.lenses[2].correction(&shot(8.0)).is_none());
|
||||
assert!(db.lenses[0].correction(&Shot { width: f64::NAN, ..shot(24.0) }).is_none());
|
||||
assert!(db.lenses[0].correction(&Shot { focal: 0.0, ..shot(24.0) }).is_none());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn load_dir_reads_xml_files_and_reports_bad_ones() {
|
||||
let dir = std::env::temp_dir().join(format!("lensdb-test-{}", std::process::id()));
|
||||
std::fs::create_dir_all(&dir).unwrap();
|
||||
std::fs::write(dir.join("a.xml"), DB).unwrap();
|
||||
std::fs::write(dir.join("b.xml"), "<lensdatabase").unwrap();
|
||||
std::fs::write(dir.join("c.txt"), "ignored").unwrap();
|
||||
let (db, errs) = LensDb::load_dir(&dir).unwrap();
|
||||
std::fs::remove_dir_all(&dir).unwrap();
|
||||
assert_eq!(db.lenses.len(), 3);
|
||||
assert_eq!(errs.len(), 1, "{errs:?}");
|
||||
assert!(LensDb::load_dir(Path::new("/nonexistent/lensdb")).is_err());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tokens_split_names_like_exif_writes_them() {
|
||||
assert_eq!(tokens("FE 24-105mm F4.0 G OSS"), ["fe", "24", "105", "mm", "f", "4", "g", "oss"]);
|
||||
assert_eq!(tokens("f/2.8"), ["f", "2.8"]);
|
||||
}
|
||||
@@ -0,0 +1,136 @@
|
||||
//! A small, bounded XML tree reader for lensfun-format database files.
|
||||
|
||||
use quick_xml::events::Event;
|
||||
|
||||
use crate::Error;
|
||||
|
||||
const MAX_NODES: usize = 500_000;
|
||||
const MAX_DEPTH: usize = 32;
|
||||
|
||||
#[derive(Debug, Default)]
|
||||
pub(crate) struct Node {
|
||||
pub name: String,
|
||||
pub attrs: Vec<(String, String)>,
|
||||
pub text: String,
|
||||
pub children: Vec<Node>,
|
||||
}
|
||||
|
||||
impl Node {
|
||||
pub fn attr(&self, k: &str) -> Option<&str> {
|
||||
self.attrs.iter().find(|(n, _)| n == k).map(|(_, v)| v.as_str())
|
||||
}
|
||||
pub fn num(&self, k: &str) -> Option<f64> {
|
||||
self.attr(k).and_then(|v| v.trim().parse::<f64>().ok()).filter(|v| v.is_finite())
|
||||
}
|
||||
pub fn kids<'a>(&'a self, name: &'a str) -> impl Iterator<Item = &'a Node> + 'a {
|
||||
self.children.iter().filter(move |c| c.name == name)
|
||||
}
|
||||
pub fn kid(&self, name: &str) -> Option<&Node> {
|
||||
self.children.iter().find(|c| c.name == name)
|
||||
}
|
||||
pub fn kid_text(&self, name: &str) -> Option<&str> {
|
||||
self.kid(name).map(|n| n.text.trim()).filter(|t| !t.is_empty())
|
||||
}
|
||||
}
|
||||
|
||||
fn entity(name: &str) -> Option<char> {
|
||||
Some(match name {
|
||||
"amp" => '&',
|
||||
"lt" => '<',
|
||||
"gt" => '>',
|
||||
"quot" => '"',
|
||||
"apos" => '\'',
|
||||
_ => {
|
||||
let n = name.strip_prefix('#')?;
|
||||
let v = match n.strip_prefix('x') {
|
||||
Some(h) => u32::from_str_radix(h, 16).ok()?,
|
||||
None => n.parse().ok()?,
|
||||
};
|
||||
char::from_u32(v)?
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
fn unescape(s: &str) -> String {
|
||||
let mut out = String::with_capacity(s.len());
|
||||
let mut rest = s;
|
||||
while let Some(i) = rest.find('&') {
|
||||
out.push_str(rest.get(..i).unwrap_or_default());
|
||||
let tail = rest.get(i..).unwrap_or_default();
|
||||
match tail.find(';').and_then(|j| entity(tail.get(1..j)?).map(|c| (c, j))) {
|
||||
Some((c, j)) => {
|
||||
out.push(c);
|
||||
rest = tail.get(j + 1..).unwrap_or_default();
|
||||
}
|
||||
None => {
|
||||
out.push('&');
|
||||
rest = tail.get(1..).unwrap_or_default();
|
||||
}
|
||||
}
|
||||
}
|
||||
out.push_str(rest);
|
||||
out
|
||||
}
|
||||
|
||||
pub(crate) fn parse(s: &str) -> Result<Node, Error> {
|
||||
let mut reader = quick_xml::Reader::from_str(s);
|
||||
reader.config_mut().check_end_names = false;
|
||||
let mut stack = vec![Node { name: "#root".into(), ..Default::default() }];
|
||||
let mut count = 0usize;
|
||||
loop {
|
||||
let ev = reader.read_event().map_err(|e| Error::Xml(e.to_string()))?;
|
||||
match ev {
|
||||
Event::Start(ref e) | Event::Empty(ref e) => {
|
||||
count += 1;
|
||||
if count > MAX_NODES || stack.len() > MAX_DEPTH {
|
||||
return Err(Error::TooLarge);
|
||||
}
|
||||
let attrs = e
|
||||
.attributes()
|
||||
.with_checks(false)
|
||||
.flatten()
|
||||
.map(|a| (String::from_utf8_lossy(a.key.as_ref()).into_owned(), unescape(&String::from_utf8_lossy(&a.value))))
|
||||
.collect();
|
||||
let node = Node { name: String::from_utf8_lossy(e.name().as_ref()).into_owned(), attrs, ..Default::default() };
|
||||
if matches!(ev, Event::Empty(_)) {
|
||||
if let Some(p) = stack.last_mut() {
|
||||
p.children.push(node);
|
||||
}
|
||||
} else {
|
||||
stack.push(node);
|
||||
}
|
||||
}
|
||||
Event::End(_) => close(&mut stack),
|
||||
Event::Text(t) => {
|
||||
if let Some(top) = stack.last_mut() {
|
||||
top.text.push_str(&t.decode().map_err(|e| Error::Xml(e.to_string()))?);
|
||||
}
|
||||
}
|
||||
Event::CData(t) => {
|
||||
if let Some(top) = stack.last_mut() {
|
||||
top.text.push_str(&String::from_utf8_lossy(&t));
|
||||
}
|
||||
}
|
||||
Event::GeneralRef(r) => {
|
||||
if let (Some(top), Some(c)) = (stack.last_mut(), entity(&String::from_utf8_lossy(&r))) {
|
||||
top.text.push(c);
|
||||
}
|
||||
}
|
||||
Event::Eof => break,
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
while stack.len() > 1 {
|
||||
close(&mut stack);
|
||||
}
|
||||
stack.pop().ok_or(Error::Xml("empty document".into()))
|
||||
}
|
||||
|
||||
fn close(stack: &mut Vec<Node>) {
|
||||
if stack.len() > 1
|
||||
&& let Some(n) = stack.pop()
|
||||
&& let Some(p) = stack.last_mut()
|
||||
{
|
||||
p.children.push(n);
|
||||
}
|
||||
}
|
||||
@@ -46,6 +46,7 @@ pub const TABLE: &[(&str, Class)] = &[
|
||||
("codecs", Class::Layer(1)),
|
||||
("meta", Class::Layer(1)),
|
||||
("develop", Class::Layer(1)),
|
||||
("lensdb", Class::Layer(2)),
|
||||
("faces", Class::Layer(1)),
|
||||
("scenes", Class::Layer(1)),
|
||||
("pipeline", Class::Layer(2)),
|
||||
|
||||
Reference in new issue
Block a user