1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
//! Higher-order distortion factors via Dual2 Hessians (T7.3).
//!
//! Provides curvature information beyond what [`crate::FactorsExact`] (first-order) gives,
//! by using [`crate::autodiff::Dual2`] to compute diagonal second partial derivatives
//! in a single forward pass through the projection.
use crate::autodiff::Dual2;
use crate::projection_generic::ProjectGeneric;
/// Second-order distortion factors computed via `Dual2<2>`.
///
/// Extends [`crate::FactorsExact`] with diagonal Hessian entries,
/// enabling curvature and non-linearity analysis of the map projection.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct HessianFactors {
/// Geodetic longitude λ in radians (input point).
pub lambda: f64,
/// Geodetic latitude φ in radians (input point).
pub phi: f64,
// First-order Jacobian entries (same as FactorsExact partial derivatives).
/// ∂x/∂λ
pub dxdlam: f64,
/// ∂x/∂φ
pub dxdphi: f64,
/// ∂y/∂λ
pub dydlam: f64,
/// ∂y/∂φ
pub dydphi: f64,
// Second-order diagonal Hessian entries.
/// ∂²x/∂λ²
pub d2xdlam2: f64,
/// ∂²x/∂φ²
pub d2xdphi2: f64,
/// ∂²y/∂λ²
pub d2ydlam2: f64,
/// ∂²y/∂φ²
pub d2ydphi2: f64,
/// Approximate curvature indicator: `|∂²x/∂λ² · ∂²y/∂φ² − ∂²x/∂φ² · ∂²y/∂λ²|`.
/// Non-zero when the mapping is non-linear; zero for purely affine maps.
pub curvature_indicator: f64,
}
impl HessianFactors {
/// Compute second-order distortion factors for a projection at (lambda, phi).
///
/// Uses `Dual2<2>` for exact first and second partial derivatives in one forward pass.
/// Returns `None` if the projection returns an error at this point.
pub fn compute<P: ProjectGeneric>(proj: &P, lambda: f64, phi: f64) -> Option<Self> {
let lam_d = Dual2::<2>::variable(lambda, 0);
let phi_d = Dual2::<2>::variable(phi, 1);
let (xd, yd) = proj.project_fwd_generic(lam_d, phi_d).ok()?;
let dxdlam = xd.d1[0];
let dxdphi = xd.d1[1];
let dydlam = yd.d1[0];
let dydphi = yd.d1[1];
let d2xdlam2 = xd.d2[0];
let d2xdphi2 = xd.d2[1];
let d2ydlam2 = yd.d2[0];
let d2ydphi2 = yd.d2[1];
let curvature_indicator = (d2xdlam2 * d2ydphi2 - d2xdphi2 * d2ydlam2).abs();
Some(Self {
lambda,
phi,
dxdlam,
dxdphi,
dydlam,
dydphi,
d2xdlam2,
d2xdphi2,
d2ydlam2,
d2ydphi2,
curvature_indicator,
})
}
/// Returns `true` when all second-order partial derivatives are smaller than `tolerance`.
///
/// An affine (locally-linear) map has zero second derivatives, so this
/// checks how far the projection departs from local linearity.
pub fn is_locally_linear(&self, tolerance: f64) -> bool {
self.d2xdlam2.abs() < tolerance
&& self.d2xdphi2.abs() < tolerance
&& self.d2ydlam2.abs() < tolerance
&& self.d2ydphi2.abs() < tolerance
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::scalar::Scalar;
// A simple stub projection: Mercator on unit sphere.
// x = lam, y = ln(tan(pi/4 + phi/2))
struct UnitMercator;
impl ProjectGeneric for UnitMercator {
fn project_fwd_generic<S: Scalar>(
&self,
lam: S,
phi: S,
) -> crate::error::ProjResult<(S, S)> {
let x = lam;
let y = (S::quarter_pi() + phi * S::from_f64(0.5)).tan().ln();
Ok((x, y))
}
fn project_inv_generic<S: Scalar>(&self, x: S, y: S) -> crate::error::ProjResult<(S, S)> {
let lam = x;
let phi = S::from_f64(2.0) * y.exp().atan() - S::half_pi();
Ok((lam, phi))
}
}
#[test]
fn dual2_second_derivative_of_x_squared() {
// f(x) = x², f'(x) = 2x, f''(x) = 2
// Using Dual2 directly (not via HessianFactors).
let x = Dual2::<1>::variable(3.0, 0);
let x2 = x * x;
assert!((x2.v - 9.0).abs() < 1e-10, "val = {}", x2.v);
assert!((x2.d1[0] - 6.0).abs() < 1e-10, "grad = {}", x2.d1[0]);
assert!((x2.d2[0] - 2.0).abs() < 1e-10, "hess_diag = {}", x2.d2[0]);
}
#[test]
fn dual2_second_derivative_of_sin() {
// f(x) = sin(x), f''(x) = -sin(x)
use core::f64::consts::PI;
let x = Dual2::<1>::variable(PI / 6.0, 0);
let fx = x.sin();
assert!((fx.v - 0.5).abs() < 1e-10, "val = {}", fx.v);
assert!(
(fx.d1[0] - (PI / 6.0).cos()).abs() < 1e-10,
"grad = {}",
fx.d1[0]
);
assert!(
(fx.d2[0] - (-(PI / 6.0).sin())).abs() < 1e-10,
"hess_diag = {}",
fx.d2[0]
);
}
#[test]
fn hessian_factors_mercator_at_equator() {
// Mercator at equator (phi=0, lam=0):
// x = lam => dx/dlam=1, dx/dphi=0, d2x/dlam2=0, d2x/dphi2=0
// y = ln(tan(pi/4 + phi/2)):
// dy/dphi = sec(phi) = 1 at phi=0
// d2y/dphi2 = sec(phi)*tan(phi) = 0 at phi=0
// curvature_indicator = |0*0 - 0*0| = 0
let proj = UnitMercator;
let hf =
HessianFactors::compute(&proj, 0.0, 0.0).expect("Mercator should not fail at equator");
assert!((hf.dxdlam - 1.0).abs() < 1e-10, "dxdlam = {}", hf.dxdlam);
assert!(hf.dxdphi.abs() < 1e-10, "dxdphi = {}", hf.dxdphi);
assert!(hf.dydlam.abs() < 1e-10, "dydlam = {}", hf.dydlam);
assert!((hf.dydphi - 1.0).abs() < 1e-10, "dydphi = {}", hf.dydphi);
assert!(hf.d2xdlam2.abs() < 1e-10, "d2xdlam2 = {}", hf.d2xdlam2);
assert!(hf.d2xdphi2.abs() < 1e-10, "d2xdphi2 = {}", hf.d2xdphi2);
assert!(hf.d2ydlam2.abs() < 1e-10, "d2ydlam2 = {}", hf.d2ydlam2);
assert!(hf.d2ydphi2.abs() < 1e-10, "d2ydphi2 = {}", hf.d2ydphi2);
assert!(
hf.curvature_indicator < 1e-10,
"K = {}",
hf.curvature_indicator
);
assert!(hf.is_locally_linear(1e-9));
}
#[test]
fn hessian_factors_mercator_at_phi30() {
// At phi=30deg=pi/6:
// dy/dphi = sec(phi), d2y/dphi2 = sec(phi)*tan(phi) (nonzero)
// d2x/... = 0 (x = lam is linear)
// curvature_indicator = |0 * nonzero - 0 * 0| = 0
// But is_locally_linear should be false (d2y/dphi2 != 0)
use core::f64::consts::FRAC_PI_6;
let phi = FRAC_PI_6;
let proj = UnitMercator;
let hf =
HessianFactors::compute(&proj, 0.0, phi).expect("Mercator should not fail at phi=30");
let sec_phi = 1.0 / phi.cos();
let expected_d2ydphi2 = sec_phi * phi.tan();
assert!((hf.dydphi - sec_phi).abs() < 1e-9, "dydphi = {}", hf.dydphi);
assert!(
(hf.d2ydphi2 - expected_d2ydphi2).abs() < 1e-9,
"d2ydphi2 = {}",
hf.d2ydphi2
);
assert!(
!hf.is_locally_linear(1e-9),
"Mercator at phi=30 is not locally linear"
);
}
#[test]
fn hessian_factors_fields_accessible() {
let proj = UnitMercator;
let hf = HessianFactors::compute(&proj, 0.1, 0.2).unwrap();
// All public fields must be accessible (compilation test).
let _ = (
hf.lambda, hf.phi, hf.dxdlam, hf.dxdphi, hf.dydlam, hf.dydphi,
);
let _ = (hf.d2xdlam2, hf.d2xdphi2, hf.d2ydlam2, hf.d2ydphi2);
let _ = hf.curvature_indicator;
}
}