rigidity_io/
las_format.rs

1//! Reading LAS and LAZ through the `las` crate.
2
3use std::path::Path;
4
5use nalgebra::Vector3;
6use rigidity_core::{Attribute, AttributeData, PointCloud};
7
8use crate::IoError;
9
10/// Reads a cloud from LAS/LAZ.
11///
12/// # Origin
13///
14/// LAS is almost always georeferenced: coordinates are UTM or a similar
15/// projection with values around a million metres. The origin is placed at
16/// the centre of the bounding box, since otherwise `f32` storage would
17/// destroy millimetres before the first computation.
18///
19/// # Attributes
20///
21/// Intensity is always carried over; colour is carried when the point
22/// format has it.
23pub fn read_las(path: &Path) -> Result<PointCloud, IoError> {
24    let mut reader = las::Reader::from_path(path).map_err(|e| IoError::Las(e.to_string()))?;
25    let data = reader.read_all().map_err(|e| IoError::Las(e.to_string()))?;
26
27    // First pass: the bounds, to choose the origin.
28    let mut min = Vector3::repeat(f64::INFINITY);
29    let mut max = Vector3::repeat(f64::NEG_INFINITY);
30    let mut count = 0usize;
31    for point in data.points() {
32        let point = point.map_err(|e| IoError::Las(e.to_string()))?;
33        let p = Vector3::new(point.x, point.y, point.z);
34        min = min.inf(&p);
35        max = max.sup(&p);
36        count += 1;
37    }
38    if count == 0 {
39        return Ok(PointCloud::new());
40    }
41    let origin = (min + max) * 0.5;
42
43    // Second pass: the points themselves and their attributes.
44    let mut cloud = PointCloud::with_origin(origin);
45    let mut intensity = Vec::with_capacity(count);
46    let mut colors: Option<(Vec<u16>, Vec<u16>, Vec<u16>)> = None;
47    for point in data.points() {
48        let point = point.map_err(|e| IoError::Las(e.to_string()))?;
49        cloud.push(Vector3::new(point.x, point.y, point.z));
50        intensity.push(point.intensity);
51        if let Some(color) = point.color {
52            let channels = colors.get_or_insert_with(|| {
53                (
54                    Vec::with_capacity(count),
55                    Vec::with_capacity(count),
56                    Vec::with_capacity(count),
57                )
58            });
59            channels.0.push(color.red);
60            channels.1.push(color.green);
61            channels.2.push(color.blue);
62        }
63    }
64
65    cloud.push_attribute(Attribute {
66        name: "intensity".into(),
67        data: AttributeData::U16(intensity),
68    })?;
69    if let Some((red, green, blue)) = colors {
70        for (name, channel) in [("red", red), ("green", green), ("blue", blue)] {
71            cloud.push_attribute(Attribute {
72                name: name.into(),
73                data: AttributeData::U16(channel),
74            })?;
75        }
76    }
77    Ok(cloud)
78}
79
80/// Writes a cloud as LAS, or as LAZ when the path says so.
81///
82/// The scale is a millimetre. LAS stores coordinates as scaled integers
83/// about an offset, and the scale is the quantum: at a tenth of a
84/// millimetre a survey-sized extent overflows the thirty-two bits the
85/// format gives each axis, and at a centimetre the file is coarser than
86/// the instrument that made it. A millimetre reaches ±2000 km from the
87/// offset, which is more than any projected coordinate system needs.
88pub fn write_las(cloud: &PointCloud, path: &Path) -> Result<(), IoError> {
89    /// Metres per stored unit.
90    const SCALE: f64 = 0.001;
91
92    let origin = cloud.origin();
93    let mut header = las::Builder::from((1, 4));
94    header.transforms = las::Vector {
95        x: las::Transform {
96            scale: SCALE,
97            offset: origin.x,
98        },
99        y: las::Transform {
100            scale: SCALE,
101            offset: origin.y,
102        },
103        z: las::Transform {
104            scale: SCALE,
105            offset: origin.z,
106        },
107    };
108    // Compression is chosen by the extension, the same way the reader
109    // chooses it.
110    header.point_format.is_compressed = path
111        .extension()
112        .is_some_and(|extension| extension.eq_ignore_ascii_case("laz"));
113    let header = header
114        .into_header()
115        .map_err(|e| IoError::Las(e.to_string()))?;
116
117    let mut writer =
118        las::Writer::from_path(path, header).map_err(|e| IoError::Las(e.to_string()))?;
119    let intensity = cloud.attribute("intensity");
120    for index in 0..cloud.len() {
121        let point = cloud.point(index);
122        writer
123            .write_point(las::Point {
124                x: point.x,
125                y: point.y,
126                z: point.z,
127                intensity: match intensity.map(|attribute| &attribute.data) {
128                    Some(AttributeData::U16(values)) => values.get(index).copied().unwrap_or(0),
129                    _ => 0,
130                },
131                ..Default::default()
132            })
133            .map_err(|e| IoError::Las(e.to_string()))?;
134    }
135    writer.close().map_err(|e| IoError::Las(e.to_string()))?;
136    Ok(())
137}
rigidity_io/las_format.rs

rigidity_io/
las_format.rs
