1mod error;
24
25use core::fmt;
26use core::ops::Mul;
27
28pub use error::AffineError;
29
30pub const DEFAULT_EPSILON: f64 = 1.0e-5;
32const DEFAULT_EPSILON_SQUARED: f64 = DEFAULT_EPSILON * DEFAULT_EPSILON;
33
34#[derive(Clone, Copy, Debug, PartialEq)]
43pub struct Affine {
44 pub a: f64,
45 pub b: f64,
46 pub c: f64,
47 pub d: f64,
48 pub e: f64,
49 pub f: f64,
50}
51
52impl Affine {
53 pub const IDENTITY: Self = Self::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0);
55
56 #[must_use]
58 #[allow(clippy::many_single_char_names)]
59 pub const fn new(a: f64, b: f64, c: f64, d: f64, e: f64, f: f64) -> Self {
60 Self { a, b, c, d, e, f }
61 }
62
63 #[must_use]
65 pub const fn translation(x_offset: f64, y_offset: f64) -> Self {
66 Self::new(1.0, 0.0, x_offset, 0.0, 1.0, y_offset)
67 }
68
69 #[must_use]
71 pub const fn scale(x_scale: f64, y_scale: f64) -> Self {
72 Self::new(x_scale, 0.0, 0.0, 0.0, y_scale, 0.0)
73 }
74
75 #[must_use]
77 pub const fn uniform_scale(scale: f64) -> Self {
78 Self::scale(scale, scale)
79 }
80
81 #[must_use]
83 pub fn shear(x_angle_degrees: f64, y_angle_degrees: f64) -> Self {
84 let x_shear = x_angle_degrees.to_radians().tan();
85 let y_shear = y_angle_degrees.to_radians().tan();
86 Self::new(1.0, x_shear, 0.0, y_shear, 1.0, 0.0)
87 }
88
89 #[must_use]
91 pub fn rotation(angle_degrees: f64) -> Self {
92 let (cosine, sine) = cos_sin_degrees(angle_degrees);
93 Self::new(cosine, -sine, 0.0, sine, cosine, 0.0)
94 }
95
96 #[must_use]
98 pub fn rotation_around(angle_degrees: f64, pivot: [f64; 2]) -> Self {
99 let (cosine, sine) = cos_sin_degrees(angle_degrees);
100 let [pivot_x, pivot_y] = pivot;
101 Self::new(
102 cosine,
103 -sine,
104 pivot_x - pivot_x * cosine + pivot_y * sine,
105 sine,
106 cosine,
107 pivot_y - pivot_x * sine - pivot_y * cosine,
108 )
109 }
110
111 #[must_use]
113 pub const fn permutation() -> Self {
114 Self::new(0.0, 1.0, 0.0, 1.0, 0.0, 0.0)
115 }
116
117 #[must_use]
119 #[allow(clippy::many_single_char_names)]
120 pub const fn from_gdal(coefficients: [f64; 6]) -> Self {
121 let [c, a, b, f, d, e] = coefficients;
122 Self::new(a, b, c, d, e, f)
123 }
124
125 #[must_use]
127 pub const fn to_gdal(self) -> [f64; 6] {
128 [self.c, self.a, self.b, self.f, self.d, self.e]
129 }
130
131 #[must_use]
133 pub const fn to_shapely(self) -> [f64; 6] {
134 [self.a, self.b, self.d, self.e, self.c, self.f]
135 }
136
137 #[must_use]
139 pub fn determinant(self) -> f64 {
140 self.a.mul_add(self.e, -(self.b * self.d))
141 }
142
143 #[must_use]
145 pub fn is_finite(self) -> bool {
146 [self.a, self.b, self.c, self.d, self.e, self.f]
147 .into_iter()
148 .all(f64::is_finite)
149 }
150
151 #[must_use]
153 pub fn is_degenerate(self) -> bool {
154 self.determinant() == 0.0
155 }
156
157 #[must_use]
159 pub fn is_proper(self) -> bool {
160 self.determinant() > 0.0
161 }
162
163 #[must_use]
165 pub const fn column_vectors(self) -> [[f64; 2]; 3] {
166 [[self.a, self.d], [self.b, self.e], [self.c, self.f]]
167 }
168
169 #[must_use]
171 pub fn is_identity(self) -> bool {
172 self.approx_eq(Self::IDENTITY, DEFAULT_EPSILON)
173 }
174
175 #[must_use]
177 pub fn is_rectilinear(self) -> bool {
178 (self.a.abs() < DEFAULT_EPSILON && self.e.abs() < DEFAULT_EPSILON)
179 || (self.d.abs() < DEFAULT_EPSILON && self.b.abs() < DEFAULT_EPSILON)
180 }
181
182 #[must_use]
184 pub fn is_conformal(self) -> bool {
185 self.a.mul_add(self.b, self.d * self.e).abs() < DEFAULT_EPSILON
186 }
187
188 #[must_use]
190 pub fn is_orthonormal(self) -> bool {
191 self.is_conformal()
192 && (1.0 - self.a.mul_add(self.a, self.d * self.d)).abs() < DEFAULT_EPSILON
193 && (1.0 - self.b.mul_add(self.b, self.e * self.e)).abs() < DEFAULT_EPSILON
194 }
195
196 #[must_use]
198 pub fn scaling(self) -> [f64; 2] {
199 let trace = self.a * self.a + self.b * self.b + self.d * self.d + self.e * self.e;
200 let determinant_squared = self.determinant().powi(2);
201 let mut delta = trace.mul_add(trace / 4.0, -determinant_squared);
202 if delta < DEFAULT_EPSILON_SQUARED {
203 delta = 0.0;
204 }
205 let root = delta.sqrt();
206 [(trace / 2.0 + root).sqrt(), (trace / 2.0 - root).sqrt()]
207 }
208
209 #[must_use]
211 pub fn eccentricity(self) -> f64 {
212 let [major, minor] = self.scaling();
213 major.mul_add(major, -(minor * minor)).sqrt() / major
214 }
215
216 pub fn rotation_angle(self) -> Result<f64, AffineError> {
223 if !self.is_finite() {
224 return Err(AffineError::NonFiniteTransform);
225 }
226 if !self.is_proper() || self.is_degenerate() {
227 return Err(AffineError::UndefinedRotation);
228 }
229 let [major, _] = self.scaling();
230 Ok((self.d / major).atan2(self.a / major).to_degrees())
231 }
232
233 #[must_use]
235 pub fn transform_point(self, [x, y]: [f64; 2]) -> [f64; 2] {
236 [
237 x.mul_add(self.a, y.mul_add(self.b, self.c)),
238 x.mul_add(self.d, y.mul_add(self.e, self.f)),
239 ]
240 }
241
242 pub fn transform_points_in_place(self, points: &mut [[f64; 2]]) {
244 for point in points {
245 *point = self.transform_point(*point);
246 }
247 }
248
249 #[must_use]
253 pub fn compose(self, rhs: Self) -> Self {
254 Self::new(
255 self.a.mul_add(rhs.a, self.b * rhs.d),
256 self.a.mul_add(rhs.b, self.b * rhs.e),
257 self.a.mul_add(rhs.c, self.b.mul_add(rhs.f, self.c)),
258 self.d.mul_add(rhs.a, self.e * rhs.d),
259 self.d.mul_add(rhs.b, self.e * rhs.e),
260 self.d.mul_add(rhs.c, self.e.mul_add(rhs.f, self.f)),
261 )
262 }
263
264 pub fn inverse(self) -> Result<Self, AffineError> {
272 if !self.is_finite() {
273 return Err(AffineError::NonFiniteTransform);
274 }
275
276 let determinant = self.determinant();
277 if determinant == 0.0 {
278 return Err(AffineError::NonInvertibleTransform);
279 }
280
281 let inverse_determinant = determinant.recip();
282 let a = self.e * inverse_determinant;
283 let b = -self.b * inverse_determinant;
284 let d = -self.d * inverse_determinant;
285 let e = self.a * inverse_determinant;
286
287 Ok(Self::new(
288 a,
289 b,
290 -self.c.mul_add(a, self.f * b),
291 d,
292 e,
293 -self.c.mul_add(d, self.f * e),
294 ))
295 }
296
297 #[must_use]
299 pub fn approx_eq(self, other: Self, epsilon: f64) -> bool {
300 let lhs = [self.a, self.b, self.c, self.d, self.e, self.f];
301 let rhs = [other.a, other.b, other.c, other.d, other.e, other.f];
302 lhs.into_iter()
303 .zip(rhs)
304 .all(|(left, right)| (left - right).abs() < epsilon)
305 }
306
307 #[allow(clippy::many_single_char_names)]
318 pub fn from_world_file(text: &str) -> Result<Self, AffineError> {
319 let values = text
320 .split_whitespace()
321 .map(str::parse::<f64>)
322 .collect::<Result<Vec<_>, _>>()
323 .map_err(|_| AffineError::InvalidWorldFile)?;
324 let [a, d, b, e, c, f] =
325 <[f64; 6]>::try_from(values).map_err(|_| AffineError::InvalidWorldFile)?;
326 let center = Self::new(a, b, c, d, e, f);
327 if !center.is_finite() {
328 return Err(AffineError::InvalidWorldFile);
329 }
330 Ok(center.compose(Self::translation(-0.5, -0.5)))
331 }
332
333 #[must_use]
335 pub fn to_world_file(self) -> String {
336 let center = self.compose(Self::translation(0.5, 0.5));
337 [center.a, center.d, center.b, center.e, center.c, center.f]
338 .into_iter()
339 .map(|value| format!("{value:?}"))
340 .collect::<Vec<_>>()
341 .join("\n")
342 + "\n"
343 }
344}
345
346impl Default for Affine {
347 fn default() -> Self {
348 Self::IDENTITY
349 }
350}
351
352impl Mul for Affine {
353 type Output = Self;
354
355 fn mul(self, rhs: Self) -> Self::Output {
356 self.compose(rhs)
357 }
358}
359
360impl Mul<[f64; 2]> for Affine {
361 type Output = [f64; 2];
362
363 fn mul(self, rhs: [f64; 2]) -> Self::Output {
364 self.transform_point(rhs)
365 }
366}
367
368impl fmt::Display for Affine {
369 fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
370 write!(
371 formatter,
372 "|{:.2},{:.2},{:.2}|\n|{:.2},{:.2},{:.2}|\n|0.00,0.00,1.00|",
373 self.a, self.b, self.c, self.d, self.e, self.f
374 )
375 }
376}
377
378fn cos_sin_degrees(angle: f64) -> (f64, f64) {
379 let normalized = angle.rem_euclid(360.0);
380 match normalized {
381 90.0 => (0.0, 1.0),
382 180.0 => (-1.0, 0.0),
383 270.0 => (0.0, -1.0),
384 _ => {
385 let (sine, cosine) = normalized.to_radians().sin_cos();
386 (cosine, sine)
387 }
388 }
389}