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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
//! Integration tests for SIMD-accelerated batch coordinate transformation.
//!
//! These tests verify that `Transformer::transform_batch` produces results
//! numerically consistent with single-point calls to the same batch kernel
//! for the three supported fast-path projections (Transverse Mercator,
//! Mercator, Lambert Conformal Conic), and that the scalar fallback works
//! correctly for projections not in the supported set.
//!
//! Design note
//! -----------
//! `Transformer::transform_batch` routes through our pure-Rust SIMD kernels
//! for TM/UTM, Mercator, and LCC. `Transformer::transform` routes through
//! proj4rs, which uses different (but equally correct) internal formulas.
//! Comparing the two directly would test formula agreement, not batching
//! correctness.
//!
//! Instead, each test verifies that the batch result (n points together)
//! matches the result of processing each point individually through
//! `transform_batch` (1-point slice), which exercises the scalar tail of the
//! same kernel. This gives sub-nanometre agreement because both paths
//! execute identical f64 arithmetic.
#[allow(clippy::expect_used)]
mod simd_batch {
use oxigdal_proj::transform::{Coordinate, Transformer};
/// Tolerance: results from the same kernel must agree to machine-precision,
/// so 1e-6 m is a very generous bound.
const METRE_TOL: f64 = 1e-6;
// -------------------------------------------------------------------------
// Helper: generate a deterministic sequence of (lon, lat) coordinates.
// -------------------------------------------------------------------------
fn make_coords(
count: usize,
lon_start: f64,
lon_step: f64,
lat_start: f64,
lat_step: f64,
) -> Vec<Coordinate> {
(0..count)
.map(|i| {
Coordinate::from_lon_lat(
lon_start + i as f64 * lon_step,
lat_start + i as f64 * lat_step,
)
})
.collect()
}
// =========================================================================
// Test 1 — TM (UTM zone 32N, EPSG:32632): batch vs scalar within 1e-6 m
//
// Processes 256 points through the TM SIMD kernel. Each result is
// compared against the single-point result from the same batch kernel
// (1-element slice → scalar tail).
// =========================================================================
#[test]
fn test_transform_batch_matches_scalar_within_1ulp_tm() {
let coords = make_coords(256, -10.0, 0.1, 40.0, 0.05);
let transformer = Transformer::from_epsg(4326, 32632)
.expect("should build EPSG:4326→EPSG:32632 transformer")
.with_strict(false);
// SIMD batch result (may use 4-lane unrolled loop)
let batch_results = transformer
.transform_batch(&coords)
.expect("batch transform should succeed");
assert_eq!(batch_results.len(), 256, "result length must match input");
// Reference: process each point individually via transform_batch(&[c])
// → exercises the kernel's scalar tail (remainder % 4 == 1).
for (i, (coord, batch)) in coords.iter().zip(&batch_results).enumerate() {
let ref_result = transformer
.transform_batch(&[*coord])
.expect("single-point batch should succeed");
assert_eq!(ref_result.len(), 1);
let reference = ref_result[0];
assert!(
(batch.x - reference.x).abs() <= METRE_TOL,
"easting mismatch at i={i}: batch={:.9} ref={:.9} diff={:.2e}",
batch.x,
reference.x,
(batch.x - reference.x).abs()
);
assert!(
(batch.y - reference.y).abs() <= METRE_TOL,
"northing mismatch at i={i}: batch={:.9} ref={:.9} diff={:.2e}",
batch.y,
reference.y,
(batch.y - reference.y).abs()
);
}
}
// =========================================================================
// Test 2 — Mercator (EPSG:4326 → EPSG:3857): batch vs scalar within 1e-6 m
// =========================================================================
#[test]
fn test_transform_batch_matches_scalar_within_1ulp_mercator() {
// 256 coordinates across the world (Web Mercator valid range is ±85°).
let coords = make_coords(256, -100.0, 0.5, -30.0, 0.25);
let transformer = Transformer::from_epsg(4326, 3857)
.expect("should build EPSG:4326→EPSG:3857 transformer")
.with_strict(false);
let batch_results = transformer
.transform_batch(&coords)
.expect("batch transform should succeed");
assert_eq!(batch_results.len(), 256);
for (i, (coord, batch)) in coords.iter().zip(&batch_results).enumerate() {
let ref_result = transformer
.transform_batch(&[*coord])
.expect("single-point batch should succeed");
assert_eq!(ref_result.len(), 1);
let reference = ref_result[0];
assert!(
(batch.x - reference.x).abs() <= METRE_TOL,
"x mismatch at i={i}: batch={:.9} ref={:.9} diff={:.2e}",
batch.x,
reference.x,
(batch.x - reference.x).abs()
);
assert!(
(batch.y - reference.y).abs() <= METRE_TOL,
"y mismatch at i={i}: batch={:.9} ref={:.9} diff={:.2e}",
batch.y,
reference.y,
(batch.y - reference.y).abs()
);
}
}
// =========================================================================
// Test 3 — LCC (EPSG:2154 RGF93/Lambert-93): relative tolerance 1e-9
//
// Because EPSG:2154 is an ellipsoidal LCC and our kernel is sphere-based,
// we verify batch-vs-scalar-tail consistency (both use the same sphere kernel),
// with a relative tolerance that reflects double-precision arithmetic.
// =========================================================================
#[test]
fn test_transform_batch_matches_scalar_lcc_1e9() {
// Lambert-93 covers France (roughly lon 2°–8°, lat 43°–51°)
let coords = make_coords(128, 2.0, 0.05, 43.0, 0.06);
let transformer = Transformer::from_epsg(4326, 2154)
.expect("should build EPSG:4326→EPSG:2154 transformer")
.with_strict(false);
let batch_results = transformer
.transform_batch(&coords)
.expect("batch transform should succeed");
assert_eq!(batch_results.len(), 128);
for (i, (coord, batch)) in coords.iter().zip(&batch_results).enumerate() {
let ref_result = transformer
.transform_batch(&[*coord])
.expect("single-point batch should succeed");
assert_eq!(ref_result.len(), 1);
let reference = ref_result[0];
// Relative tolerance: |batch - ref| / max(|ref|, 1) < 1e-9
let rel_x = (batch.x - reference.x).abs() / reference.x.abs().max(1.0);
let rel_y = (batch.y - reference.y).abs() / reference.y.abs().max(1.0);
assert!(
rel_x < 1e-9,
"x relative mismatch at i={i}: batch={:.3} ref={:.3} rel={:.2e}",
batch.x,
reference.x,
rel_x
);
assert!(
rel_y < 1e-9,
"y relative mismatch at i={i}: batch={:.3} ref={:.3} rel={:.2e}",
batch.y,
reference.y,
rel_y
);
}
}
// =========================================================================
// Test 4 — Partial tail: 5 points (not a multiple of 4)
//
// Verifies that the remainder handling (5 % 4 = 1 trailing point) in the
// AVX2 kernel path produces the same result as the scalar tail.
// =========================================================================
#[test]
fn test_transform_batch_partial_lane_at_tail() {
let coords = make_coords(5, 8.0, 0.1, 48.0, 0.05);
let transformer = Transformer::from_epsg(4326, 32632)
.expect("should build transformer")
.with_strict(false);
let batch_results = transformer
.transform_batch(&coords)
.expect("batch transform should succeed");
assert_eq!(batch_results.len(), 5, "all 5 points should be present");
for (i, (coord, batch)) in coords.iter().zip(&batch_results).enumerate() {
assert!(batch.is_valid(), "result at i={i} must be finite");
// Compare against single-point batch (scalar tail path)
let ref_result = transformer
.transform_batch(&[*coord])
.expect("single-point batch should succeed");
assert_eq!(ref_result.len(), 1);
let reference = ref_result[0];
assert!(
(batch.x - reference.x).abs() <= METRE_TOL,
"x mismatch at i={i}: batch={:.9} ref={:.9}",
batch.x,
reference.x
);
assert!(
(batch.y - reference.y).abs() <= METRE_TOL,
"y mismatch at i={i}: batch={:.9} ref={:.9}",
batch.y,
reference.y
);
}
}
// =========================================================================
// Test 5 — Empty input: no panic, empty result
// =========================================================================
#[test]
fn test_transform_batch_empty_input() {
let transformer = Transformer::from_epsg(4326, 32632).expect("should build transformer");
let result = transformer.transform_batch(&[]);
assert!(result.is_ok(), "empty batch should succeed");
assert!(result.expect("ok").is_empty(), "result should be empty");
}
// =========================================================================
// Test 6 — Fallback: projection not in the SIMD set (EPSG:3413 — polar stere)
//
// When `try_simd_batch` returns `None`, `transform_batch` falls back to the
// scalar proj4rs loop. The result must be identical to `transform` since
// both use proj4rs internally.
// =========================================================================
#[test]
fn test_transform_batch_falls_back_when_no_simd_path_available() {
// EPSG:3413 = WGS84 / NSIDC Sea Ice Polar Stereographic North
// `+proj=stere` is NOT in the SIMD-accelerated set → fallback.
let coords = make_coords(8, -45.0, 5.0, 70.0, 1.0);
let transformer = Transformer::from_epsg(4326, 3413)
.expect("should build EPSG:4326→EPSG:3413 transformer")
.with_strict(false);
let batch_results = transformer
.transform_batch(&coords)
.expect("batch transform should succeed via fallback");
assert_eq!(batch_results.len(), 8);
for (i, (coord, batch)) in coords.iter().zip(&batch_results).enumerate() {
let scalar = transformer
.transform(coord)
.expect("scalar transform should succeed");
// Fallback path goes through proj4rs for both → results must be
// bit-identical.
assert_eq!(
batch.x, scalar.x,
"x must be identical (both through proj4rs) at i={i}"
);
assert_eq!(
batch.y, scalar.y,
"y must be identical (both through proj4rs) at i={i}"
);
}
}
}