pasture_algorithms/
reprojection.rs

1use std::ffi::CString;
2
3use anyhow::Result;
4use pasture_core::containers::{BorrowedBufferExt, BorrowedMutBuffer, BorrowedMutBufferExt};
5use pasture_core::math::AABB;
6use pasture_core::nalgebra::{Point3, Vector3};
7
8use pasture_core::layout::attributes::POSITION_3D;
9/// Wrapper around the proj types from the proj_sys crate. Supports transformations (the Rust proj bindings don't support this)
10pub struct Projection {
11    proj_context: *mut proj_sys::pj_ctx,
12    projection: *mut proj_sys::PJconsts,
13}
14
15impl Projection {
16    pub fn new(source_crs: &str, target_crs: &str) -> Result<Self> {
17        let src_cstr = CString::new(source_crs)?;
18        let target_cstr = CString::new(target_crs)?;
19
20        unsafe {
21            let proj_context = proj_sys::proj_context_create();
22
23            let projection = proj_sys::proj_create_crs_to_crs(
24                proj_context,
25                src_cstr.as_ptr(),
26                target_cstr.as_ptr(),
27                std::ptr::null_mut(),
28            );
29
30            Ok(Self {
31                proj_context,
32                projection,
33            })
34        }
35    }
36
37    /// Performs a transformation of the given position
38    pub fn transform(&self, position: Vector3<f64>) -> Vector3<f64> {
39        unsafe {
40            let coord = proj_sys::proj_coord(position.x, position.y, position.z, 0.0);
41            let target_coords =
42                proj_sys::proj_trans(self.projection, proj_sys::PJ_DIRECTION_PJ_FWD, coord);
43            Vector3::new(target_coords.v[0], target_coords.v[1], target_coords.v[2])
44        }
45    }
46
47    /// Transforms the given bounding box using the associated `Projection`. This keeps the bounding box axis-aligned,
48    /// the size might change though
49    pub fn transform_bounds(&self, bounds: &AABB<f64>) -> AABB<f64> {
50        let transformed_min =
51            self.transform(Vector3::new(bounds.min().x, bounds.min().y, bounds.min().z));
52        let transformed_max =
53            self.transform(Vector3::new(bounds.max().x, bounds.max().y, bounds.max().z));
54
55        let bounds: AABB<f64> = AABB::from_min_max_unchecked(
56            Point3::from(transformed_min),
57            Point3::from(transformed_min),
58        );
59        AABB::extend_with_point(&bounds, &Point3::from(transformed_max))
60    }
61}
62
63impl Drop for Projection {
64    fn drop(&mut self) {
65        unsafe {
66            proj_sys::proj_destroy(self.projection);
67            proj_sys::proj_context_destroy(self.proj_context);
68        }
69    }
70}
71
72/// Reprojection Algorithm
73/// Rewrites the 3D coordinates from the given point cloud to the given target coordinate reference system.
74/// It iterates over all points in the given point cloud.
75/// Make sure that source_crs and target_crs are valid coordinate reference systems.
76///
77/// # Panics
78///
79/// Panics if the PointLayout of this buffer does not contain the given attribute.
80///
81/// # Examples
82///
83/// ```
84/// # use pasture_algorithms::reprojection::reproject_point_cloud_within;
85/// # use pasture_core::containers::*;
86/// # use pasture_core::layout::PointType;
87/// # use pasture_core::nalgebra::Vector3;
88/// # use pasture_derive::PointType;
89
90/// #[repr(C, packed)]
91/// #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
92/// struct SimplePoint {
93///     #[pasture(BUILTIN_POSITION_3D)]
94///     pub position: Vector3<f64>,
95///     #[pasture(BUILTIN_INTENSITY)]
96///     pub intensity: u16,
97/// }
98
99/// fn main() {
100///     let points = vec![
101///         SimplePoint {
102///             position: Vector3::new(1.0, 22.0, 0.0),
103///             intensity: 42,
104///         },
105///         SimplePoint {
106///             position: Vector3::new(12.0, 23.0, 0.0),
107///             intensity: 84,
108///         },
109///         SimplePoint {
110///             position: Vector3::new(10.0, 8.0, 2.0),
111///             intensity: 84,
112///         },
113///         SimplePoint {
114///             position: Vector3::new(10.0, 0.0, 1.0),
115///             intensity: 84,
116///         },
117///     ];
118
119///     let mut interleaved = points.into_iter().collect::<VectorBuffer>();
120
121///     reproject_point_cloud_within(
122///         &mut interleaved,
123///         "EPSG:4326",
124///         "EPSG:3309",
125///     );
126
127///     for point in interleaved.view::<SimplePoint>() {
128///         println!("{:?}", point);
129///     }
130/// }
131/// ```
132pub fn reproject_point_cloud_within<'a, T: BorrowedMutBuffer<'a>>(
133    point_cloud: &'a mut T,
134    source_crs: &str,
135    target_crs: &str,
136) {
137    let proj = Projection::new(source_crs, target_crs).unwrap();
138
139    let num_points = point_cloud.len();
140    let mut positions_view = point_cloud.view_attribute_mut(&POSITION_3D);
141    for index in 0..num_points {
142        let point = positions_view.at(index);
143        let reproj = proj.transform(point);
144        positions_view.set_at(index, reproj);
145    }
146}
147
148/// Reprojection Algorithm
149/// Rewrites the 3D coordinates from the given point cloud to the given target coordinate reference system.
150/// It iterates over all points in the given point cloud.
151/// Make sure that source_crs and target_crs are valid coordinate reference systems.
152///
153/// # Panics
154///
155/// Panics if the PointLayout of this buffer does not contain the given attribute.
156///
157/// # Examples
158///
159/// ```
160/// # use pasture_algorithms::reprojection::reproject_point_cloud_between;
161/// # use pasture_core::containers::*;
162/// # use pasture_core::layout::PointType;
163/// # use pasture_core::nalgebra::Vector3;
164/// # use pasture_derive::PointType;
165/// #[repr(C, packed)]
166/// #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
167/// struct SimplePoint {
168///     #[pasture(BUILTIN_POSITION_3D)]
169///     pub position: Vector3<f64>,
170///     #[pasture(BUILTIN_INTENSITY)]
171///     pub intensity: u16,
172/// }
173/// fn main() {
174///     let points = vec![
175///         SimplePoint {
176///             position: Vector3::new(1.0, 22.0, 0.0),
177///             intensity: 42,
178///         },
179///         SimplePoint {
180///             position: Vector3::new(12.0, 23.0, 0.0),
181///             intensity: 84,
182///         },
183///         SimplePoint {
184///             position: Vector3::new(10.0, 8.0, 2.0),
185///             intensity: 84,
186///         },
187///         SimplePoint {
188///             position: Vector3::new(10.0, 0.0, 1.0),
189///             intensity: 84,
190///         },
191///     ];
192///     let mut interleaved = points.into_iter().collect::<VectorBuffer>();
193///     let mut attribute = HashMapBuffer::with_capacity(interleaved.len(), SimplePoint::layout());
194///     attribute.resize(interleaved.len());
195///     reproject_point_cloud_between(&mut interleaved, &mut attribute, "EPSG:4326", "EPSG:3309");
196///     for point in attribute.view::<SimplePoint>() {
197///         println!("{:?}", point);
198///     }
199/// }
200/// ```
201pub fn reproject_point_cloud_between<
202    'a,
203    'b,
204    T1: BorrowedMutBuffer<'a>,
205    T2: BorrowedMutBuffer<'b>,
206>(
207    source_point_cloud: &'a mut T1,
208    target_point_cloud: &'b mut T2,
209    source_crs: &str,
210    target_crs: &str,
211) {
212    if source_point_cloud.len() != target_point_cloud.len() {
213        panic!("The point clouds don't have the same size!");
214    }
215
216    let proj = Projection::new(source_crs, target_crs).unwrap();
217
218    let mut target_positions = target_point_cloud.view_attribute_mut(&POSITION_3D);
219    for (index, point) in source_point_cloud
220        .view_attribute::<Vector3<f64>>(&POSITION_3D)
221        .into_iter()
222        .enumerate()
223    {
224        let reproj = proj.transform(point);
225        target_positions.set_at(index, reproj);
226    }
227}
228
229#[cfg(test)]
230mod tests {
231    use assert_approx_eq::assert_approx_eq;
232    use pasture_core::{
233        containers::{BorrowedBuffer, HashMapBuffer, OwningBuffer, VectorBuffer},
234        layout::PointType,
235        nalgebra::Vector3,
236    };
237    use pasture_derive::PointType;
238
239    use super::*;
240
241    #[repr(C, packed)]
242    #[derive(PointType, Debug, Clone, Copy, bytemuck::AnyBitPattern, bytemuck::NoUninit)]
243    pub struct SimplePoint {
244        #[pasture(BUILTIN_POSITION_3D)]
245        pub position: Vector3<f64>,
246        #[pasture(BUILTIN_INTENSITY)]
247        pub intensity: u16,
248    }
249
250    #[test]
251    fn reproject_epsg4326_epsg3309_within() {
252        let points = vec![
253            SimplePoint {
254                position: Vector3::new(1.0, 22.0, 0.0),
255                intensity: 42,
256            },
257            SimplePoint {
258                position: Vector3::new(12.0, 23.0, 0.0),
259                intensity: 84,
260            },
261            SimplePoint {
262                position: Vector3::new(10.0, 8.0, 2.0),
263                intensity: 84,
264            },
265            SimplePoint {
266                position: Vector3::new(10.0, 0.0, 1.0),
267                intensity: 84,
268            },
269        ];
270
271        let mut interleaved = points.into_iter().collect::<VectorBuffer>();
272
273        reproject_point_cloud_within(&mut interleaved, "EPSG:4326", "EPSG:3309");
274
275        let results = [
276            Vector3::new(12185139.590523569, 7420953.944297638, 0.0),
277            Vector3::new(11104667.534080556, 7617693.973680517, 0.0),
278            Vector3::new(11055663.927418157, 5832081.512011217, 2.0),
279            Vector3::new(10807262.110686881, 4909128.916889962, 1.0),
280        ];
281
282        for (index, coord) in interleaved
283            .view_attribute::<Vector3<f64>>(&POSITION_3D)
284            .into_iter()
285            .enumerate()
286        {
287            assert_approx_eq!(coord[0], results[index][0], 0.0001);
288            assert_approx_eq!(coord[1], results[index][1], 0.0001);
289            assert_approx_eq!(coord[2], results[index][2], 0.0001);
290        }
291    }
292    #[test]
293    fn reproject_epsg4326_epsg3309_between() {
294        let points = vec![
295            SimplePoint {
296                position: Vector3::new(1.0, 22.0, 0.0),
297                intensity: 42,
298            },
299            SimplePoint {
300                position: Vector3::new(12.0, 23.0, 0.0),
301                intensity: 84,
302            },
303            SimplePoint {
304                position: Vector3::new(10.0, 8.0, 2.0),
305                intensity: 84,
306            },
307            SimplePoint {
308                position: Vector3::new(10.0, 0.0, 1.0),
309                intensity: 84,
310            },
311        ];
312
313        let mut interleaved = points.into_iter().collect::<VectorBuffer>();
314
315        let mut attribute = HashMapBuffer::with_capacity(interleaved.len(), SimplePoint::layout());
316
317        attribute.resize(interleaved.len());
318
319        reproject_point_cloud_between(&mut interleaved, &mut attribute, "EPSG:4326", "EPSG:3309");
320
321        let results = [
322            Vector3::new(12185139.590523569, 7420953.944297638, 0.0),
323            Vector3::new(11104667.534080556, 7617693.973680517, 0.0),
324            Vector3::new(11055663.927418157, 5832081.512011217, 2.0),
325            Vector3::new(10807262.110686881, 4909128.916889962, 1.0),
326        ];
327
328        for (index, coord) in attribute
329            .view_attribute::<Vector3<f64>>(&POSITION_3D)
330            .into_iter()
331            .enumerate()
332        {
333            assert_approx_eq!(coord[0], results[index][0], 0.0001);
334            assert_approx_eq!(coord[1], results[index][1], 0.0001);
335            assert_approx_eq!(coord[2], results[index][2], 0.0001);
336        }
337    }
338    #[test]
339    #[should_panic(expected = "The point clouds don't have the same size!")]
340    fn reproject_epsg4326_epsg3309_between_error() {
341        let points = vec![
342            SimplePoint {
343                position: Vector3::new(1.0, 22.0, 0.0),
344                intensity: 42,
345            },
346            SimplePoint {
347                position: Vector3::new(12.0, 23.0, 0.0),
348                intensity: 84,
349            },
350            SimplePoint {
351                position: Vector3::new(10.0, 8.0, 2.0),
352                intensity: 84,
353            },
354            SimplePoint {
355                position: Vector3::new(10.0, 0.0, 1.0),
356                intensity: 84,
357            },
358        ];
359
360        let mut interleaved = points.into_iter().collect::<VectorBuffer>();
361
362        let mut attribute = HashMapBuffer::with_capacity(2, SimplePoint::layout());
363
364        attribute.resize(2);
365
366        reproject_point_cloud_between(&mut interleaved, &mut attribute, "EPSG:4326", "EPSG:3309");
367    }
368}
pasture_algorithms/reprojection.rs

pasture_algorithms/
reprojection.rs
