ogeom_io/
cloud.rs

1//! Point clouds, and the text formats that carry them: XYZ and PTS.
2//!
3//! A cloud is points with what a scanner measured at each: a normal, a
4//! colour, and any number of named scalars (intensity, confidence, a
5//! deviation an inspection computed). [`PointCloud`] holds them side by
6//! side, one entry per point in each.
7//!
8//! XYZ is a row of numbers per point, separated by spaces, tabs, commas or
9//! semicolons; PTS is the same with a line before each scan saying how many
10//! rows follow. Neither says what its columns are, so a caller may name
11//! them ([`XyzColumn`]); otherwise a header line naming them (`# x y z nx ny
12//! nz`, `//X,Y,Z,R,G,B`) is honoured, and otherwise the column count decides,
13//! as [`read_xyz`] lays out, and a count that leaves the reading open is
14//! refused by name. Neither declares a unit; coordinates come back as
15//! written.
16
17use ogeom_core::{OgeomResult, ogeom_bail};
18use ogeom_math::{Point, Vector};
19use std::fmt::Write as _;
20
21/// One named value per point.
22#[derive(Debug, Clone, PartialEq)]
23pub struct Scalar {
24    /// The name the file gives it.
25    pub name: String,
26    /// One value per point, in point order.
27    pub values: Vec<f64>,
28}
29
30/// Points, with what was measured at each.
31///
32/// Every attribute present has exactly one entry per point;
33/// [`PointCloud::check`] says whether one does not.
34#[derive(Debug, Clone, Default, PartialEq)]
35pub struct PointCloud {
36    /// The points.
37    pub points: Vec<Point>,
38    /// A normal per point, as the file gives it, where it gives one.
39    pub normals: Option<Vec<Vector>>,
40    /// An RGBA colour per point in `[0, 1]`, where the file gives one.
41    pub colours: Option<Vec<[f64; 4]>>,
42    /// Named scalars per point, in the order the file gives them.
43    pub scalars: Vec<Scalar>,
44}
45
46impl PointCloud {
47    /// Bare points, with nothing measured at them.
48    #[must_use]
49    pub fn new(points: Vec<Point>) -> Self {
50        Self {
51            points,
52            ..Self::default()
53        }
54    }
55
56    /// The number of points.
57    #[must_use]
58    pub fn len(&self) -> usize {
59        self.points.len()
60    }
61
62    /// Whether there are no points.
63    #[must_use]
64    pub fn is_empty(&self) -> bool {
65        self.points.is_empty()
66    }
67
68    /// The values of the scalar named `name`.
69    #[must_use]
70    pub fn scalar(&self, name: &str) -> Option<&[f64]> {
71        self.scalars
72            .iter()
73            .find(|s| s.name == name)
74            .map(|s| s.values.as_slice())
75    }
76
77    /// Check that every attribute has one entry per point.
78    ///
79    /// # Errors
80    ///
81    /// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) naming
82    /// the attribute whose length differs from the point count.
83    pub fn check(&self) -> OgeomResult<()> {
84        let n = self.points.len();
85        if let Some(normals) = &self.normals
86            && normals.len() != n
87        {
88            ogeom_bail!(Dimension, "{} normals for {n} points", normals.len());
89        }
90        if let Some(colours) = &self.colours
91            && colours.len() != n
92        {
93            ogeom_bail!(Dimension, "{} colours for {n} points", colours.len());
94        }
95        for s in &self.scalars {
96            if s.values.len() != n {
97                ogeom_bail!(
98                    Dimension,
99                    "{} values of {} for {n} points",
100                    s.values.len(),
101                    s.name
102                );
103            }
104        }
105        Ok(())
106    }
107}
108
109/// What one column of an XYZ or PTS row holds.
110#[derive(Debug, Clone, PartialEq, Eq)]
111pub enum XyzColumn {
112    /// The point's `x`.
113    X,
114    /// The point's `y`.
115    Y,
116    /// The point's `z`.
117    Z,
118    /// The normal's `x`.
119    NormalX,
120    /// The normal's `y`.
121    NormalY,
122    /// The normal's `z`.
123    NormalZ,
124    /// Red, from 0 to 255.
125    Red,
126    /// Green, from 0 to 255.
127    Green,
128    /// Blue, from 0 to 255.
129    Blue,
130    /// A named scalar.
131    Scalar(String),
132    /// A column to pass over.
133    Skip,
134}
135
136impl XyzColumn {
137    /// The column a header word names; anything not a coordinate, a normal
138    /// or a colour channel is a scalar of that name.
139    fn named(word: &str) -> Self {
140        match word.to_ascii_lowercase().as_str() {
141            "x" => Self::X,
142            "y" => Self::Y,
143            "z" => Self::Z,
144            "nx" => Self::NormalX,
145            "ny" => Self::NormalY,
146            "nz" => Self::NormalZ,
147            "r" | "red" => Self::Red,
148            "g" | "green" => Self::Green,
149            "b" | "blue" => Self::Blue,
150            _ => Self::Scalar(word.to_owned()),
151        }
152    }
153
154    fn word(&self) -> &str {
155        match self {
156            Self::X => "x",
157            Self::Y => "y",
158            Self::Z => "z",
159            Self::NormalX => "nx",
160            Self::NormalY => "ny",
161            Self::NormalZ => "nz",
162            Self::Red => "red",
163            Self::Green => "green",
164            Self::Blue => "blue",
165            Self::Scalar(name) => name,
166            Self::Skip => "-",
167        }
168    }
169}
170
171fn fields(line: &str) -> impl Iterator<Item = &str> {
172    line.split(|c: char| c.is_whitespace() || c == ',' || c == ';')
173        .filter(|f| !f.is_empty())
174}
175
176/// A header line naming its columns, `x y z` first.
177fn header(line: &str) -> Option<Vec<XyzColumn>> {
178    let rest = line
179        .strip_prefix("//")
180        .or_else(|| line.strip_prefix('#'))?
181        .trim_start_matches(|c: char| c == '#' || c == '/' || c.is_whitespace());
182    let rest = rest.strip_prefix("columns:").map_or(rest, str::trim_start);
183    let columns: Vec<XyzColumn> = fields(rest).map(XyzColumn::named).collect();
184    columns
185        .starts_with(&[XyzColumn::X, XyzColumn::Y, XyzColumn::Z])
186        .then_some(columns)
187}
188
189/// The rows of a text cloud, column by column.
190struct Table {
191    columns: Vec<Vec<f64>>,
192    header: Option<Vec<XyzColumn>>,
193}
194
195/// Read the numeric rows. `counted` is PTS: a line of one integer opens a
196/// scan of that many rows.
197fn table(text: &str, counted: bool) -> OgeomResult<Table> {
198    let mut columns: Vec<Vec<f64>> = Vec::new();
199    let mut header_line = None;
200    let mut seen_data = false;
201    // PTS: the rows the open scan still owes.
202    let mut owed: Option<usize> = None;
203    for (number, line) in text.lines().enumerate() {
204        let line = line.trim();
205        let number = number + 1;
206        if line.is_empty() {
207            continue;
208        }
209        if line.starts_with('#') || line.starts_with("//") {
210            if !seen_data && header_line.is_none() {
211                header_line = header(line);
212            }
213            continue;
214        }
215        let mut row = fields(line).peekable();
216        if counted && let Some(only) = row.peek().copied() {
217            let mut probe = fields(line);
218            probe.next();
219            if probe.next().is_none() {
220                let Ok(count) = only.parse::<usize>() else {
221                    ogeom_bail!(
222                        Construction,
223                        "line {number} of the PTS holds one field, {only}, which is not a point count"
224                    );
225                };
226                if let Some(left) = owed.filter(|&n| n > 0) {
227                    ogeom_bail!(
228                        Construction,
229                        "line {number} opens a scan while the one before still owes {left} points"
230                    );
231                }
232                owed = Some(count);
233                continue;
234            }
235            match owed.as_mut() {
236                None => ogeom_bail!(
237                    Construction,
238                    "the PTS has no point count before its first point, on line {number}"
239                ),
240                Some(0) => ogeom_bail!(
241                    Construction,
242                    "line {number} is a point past the count its scan gave"
243                ),
244                Some(n) => *n -= 1,
245            }
246        }
247        let mut width = 0;
248        for (k, field) in row.enumerate() {
249            let Ok(value) = field.parse::<f64>() else {
250                ogeom_bail!(
251                    Construction,
252                    "line {number} carries {field}, which is not a number"
253                );
254            };
255            if !value.is_finite() {
256                ogeom_bail!(
257                    Construction,
258                    "line {number} carries {field}, which is not finite"
259                );
260            }
261            if !seen_data {
262                columns.push(Vec::new());
263            } else if k >= columns.len() {
264                ogeom_bail!(
265                    Construction,
266                    "line {number} has more columns than the {} the first point has",
267                    columns.len()
268                );
269            }
270            columns[k].push(value);
271            width = k + 1;
272        }
273        if seen_data && width != columns.len() {
274            ogeom_bail!(
275                Construction,
276                "line {number} has {width} columns where the first point has {}",
277                columns.len()
278            );
279        }
280        seen_data = true;
281    }
282    if let Some(left) = owed.filter(|&n| n > 0) {
283        ogeom_bail!(
284            Construction,
285            "the PTS ends with its last scan {left} points short of its count"
286        );
287    }
288    Ok(Table {
289        columns,
290        header: header_line,
291    })
292}
293
294/// Whether a column triple reads as 0-to-255 colour channels rather than a
295/// normal: every value an integer in `[0, 255]`, and one of them past the
296/// `1` no component of a unit normal reaches.
297fn colour_like(columns: &[Vec<f64>]) -> bool {
298    let channel = |v: &f64| v.fract() == 0.0 && (0.0..=255.0).contains(v);
299    columns.iter().all(|c| c.iter().all(channel))
300        && columns.iter().any(|c| c.iter().any(|&v| v > 1.0))
301}
302
303/// The columns a cloud of `width` unnamed columns is read with.
304fn inferred(columns: &[Vec<f64>]) -> OgeomResult<Vec<XyzColumn>> {
305    use XyzColumn::{Blue, Green, NormalX, NormalY, NormalZ, Red, X, Y, Z};
306    let intensity = || XyzColumn::Scalar("intensity".to_owned());
307    let colour = [Red, Green, Blue];
308    let normal = [NormalX, NormalY, NormalZ];
309    let mut layout = vec![X, Y, Z];
310    match columns.len() {
311        3 => {}
312        4 => layout.push(intensity()),
313        6 => layout.extend(if colour_like(&columns[3..6]) {
314            colour
315        } else {
316            normal
317        }),
318        7 => {
319            layout.push(intensity());
320            layout.extend(colour);
321        }
322        9 => match (colour_like(&columns[3..6]), colour_like(&columns[6..9])) {
323            (true, false) => {
324                layout.extend(colour);
325                layout.extend(normal);
326            }
327            (false, true) => {
328                layout.extend(normal);
329                layout.extend(colour);
330            }
331            _ => ogeom_bail!(
332                Construction,
333                "a cloud of nine columns reads as colour and normals only when one triple is colour; name its columns"
334            ),
335        },
336        n => ogeom_bail!(
337            Construction,
338            "a cloud of {n} columns does not say what they are; name its columns"
339        ),
340    }
341    Ok(layout)
342}
343
344fn assemble(table: Table, named: Option<&[XyzColumn]>) -> OgeomResult<PointCloud> {
345    let Table { columns, header } = table;
346    if columns.is_empty() {
347        ogeom_bail!(Construction, "the file carries no points");
348    }
349    let layout = match (named, header) {
350        (Some(named), _) => named.to_vec(),
351        (None, Some(header)) => header,
352        (None, None) => inferred(&columns)?,
353    };
354    if layout.len() != columns.len() {
355        ogeom_bail!(
356            Construction,
357            "the rows have {} columns where {} are named",
358            columns.len(),
359            layout.len()
360        );
361    }
362    let find = |wanted: &XyzColumn| -> OgeomResult<Option<usize>> {
363        let mut at = layout.iter().enumerate().filter(|(_, c)| *c == wanted);
364        match (at.next(), at.next()) {
365            (Some((i, _)), None) => Ok(Some(i)),
366            (None, _) => Ok(None),
367            (Some(_), Some(_)) => {
368                ogeom_bail!(Construction, "the column {} is named twice", wanted.word())
369            }
370        }
371    };
372    let triple = |names: [XyzColumn; 3], what: &str| -> OgeomResult<Option<[usize; 3]>> {
373        match (find(&names[0])?, find(&names[1])?, find(&names[2])?) {
374            (Some(a), Some(b), Some(c)) => Ok(Some([a, b, c])),
375            (None, None, None) => Ok(None),
376            _ => ogeom_bail!(Construction, "the columns name part of a {what}"),
377        }
378    };
379    let Some(xyz) = triple([XyzColumn::X, XyzColumn::Y, XyzColumn::Z], "point")? else {
380        ogeom_bail!(Construction, "no columns are named x, y and z");
381    };
382    let nxyz = triple(
383        [XyzColumn::NormalX, XyzColumn::NormalY, XyzColumn::NormalZ],
384        "normal",
385    )?;
386    let rgb = triple(
387        [XyzColumn::Red, XyzColumn::Green, XyzColumn::Blue],
388        "colour",
389    )?;
390    let n = columns[0].len();
391    let points = (0..n)
392        .map(|i| Point::new(columns[xyz[0]][i], columns[xyz[1]][i], columns[xyz[2]][i]))
393        .collect();
394    let normals = nxyz.map(|[a, b, c]| {
395        (0..n)
396            .map(|i| Vector::new(columns[a][i], columns[b][i], columns[c][i]))
397            .collect()
398    });
399    let colours = match rgb {
400        Some(channels) => {
401            for &c in &channels {
402                if let Some(v) = columns[c].iter().find(|v| !(0.0..=255.0).contains(*v)) {
403                    ogeom_bail!(Construction, "a colour channel reads {v}, outside 0 to 255");
404                }
405            }
406            Some(
407                (0..n)
408                    .map(|i| {
409                        let [r, g, b] = channels.map(|c| columns[c][i] / 255.0);
410                        [r, g, b, 1.0]
411                    })
412                    .collect(),
413            )
414        }
415        None => None,
416    };
417    let mut scalars: Vec<Scalar> = Vec::new();
418    let mut columns = columns;
419    for (k, column) in layout.iter().enumerate() {
420        if let XyzColumn::Scalar(name) = column {
421            if scalars.iter().any(|s| &s.name == name) {
422                ogeom_bail!(Construction, "the column {name} is named twice");
423            }
424            scalars.push(Scalar {
425                name: name.clone(),
426                values: std::mem::take(&mut columns[k]),
427            });
428        }
429    }
430    Ok(PointCloud {
431        points,
432        normals,
433        colours,
434        scalars,
435    })
436}
437
438/// Read an XYZ text cloud.
439///
440/// A row per point; blank lines and lines opening with `#` or `//` are
441/// passed over, but for the first such line before the data when it names
442/// the columns, `x y z` first: `nx ny nz` are a normal, `r g b` (or `red
443/// green blue`) a colour from 0 to 255, and any other word a scalar of
444/// that name. `columns`, when given, names them instead. Unnamed, the
445/// column count decides:
446///
447/// | columns | read as |
448/// |---------|---------|
449/// | 3 | `x y z` |
450/// | 4 | `x y z`, then a scalar `intensity` |
451/// | 6 | `x y z`, then a colour if every value of the last three is an integer from 0 to 255 and one is past 1, a normal otherwise |
452/// | 7 | `x y z`, `intensity`, then a colour, as PTS lays them out |
453/// | 9 | `x y z`, then a colour and a normal, the colour being whichever triple reads as one by the rule above |
454///
455/// Every other count is refused, as is a nine-column file where neither
456/// triple, or both, reads as colour.
457///
458/// # Errors
459///
460/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) for a
461/// field that is not a finite number, rows of differing widths, a colour
462/// channel outside 0 to 255, columns naming part of a point, normal or
463/// colour, or a layout the rule above leaves open.
464pub fn read_xyz(text: &str, columns: Option<&[XyzColumn]>) -> OgeomResult<PointCloud> {
465    assemble(table(text, false)?, columns)
466}
467
468/// Read a PTS text cloud: XYZ rows, each scan opened by a line holding its
469/// point count, the rows commonly `x y z intensity red green blue`.
470///
471/// The columns are read as [`read_xyz`] reads them; several scans in one
472/// file are read one after the other into one cloud.
473///
474/// # Errors
475///
476/// As [`read_xyz`], and for a point before any count, a scan with more or
477/// fewer rows than its count, or a count line that is not a count.
478pub fn read_pts(text: &str, columns: Option<&[XyzColumn]>) -> OgeomResult<PointCloud> {
479    assemble(table(text, true)?, columns)
480}
481
482/// Write a cloud as XYZ text.
483///
484/// A header line naming the columns, then a row per point: `x y z`, the
485/// normal, the colour as 0-to-255 integers, then each scalar, every number
486/// in the shortest form that reads back as the same `f64`. Alpha is not
487/// written; the format has no column for it. [`read_xyz`] reads the
488/// header, so the file comes back as the cloud it was, colours rounded to
489/// the nearest 255th.
490///
491/// # Errors
492///
493/// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) if an
494/// attribute does not have one entry per point (see [`PointCloud::check`]),
495/// and [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) for a scalar
496/// name that is empty, holds a separator or a `#`, or is one the header
497/// gives to a coordinate, normal or colour.
498pub fn write_xyz(cloud: &PointCloud) -> OgeomResult<String> {
499    cloud.check()?;
500    let mut out = String::from("# x y z");
501    if cloud.normals.is_some() {
502        out.push_str(" nx ny nz");
503    }
504    if cloud.colours.is_some() {
505        out.push_str(" red green blue");
506    }
507    for s in &cloud.scalars {
508        let separated = s
509            .name
510            .chars()
511            .any(|c| c.is_whitespace() || c == ',' || c == ';' || c == '#');
512        if s.name.is_empty()
513            || separated
514            || !matches!(XyzColumn::named(&s.name), XyzColumn::Scalar(_))
515        {
516            ogeom_bail!(
517                Domain,
518                "the scalar name {:?} cannot head an XYZ column",
519                s.name
520            );
521        }
522        let _ = write!(out, " {}", s.name);
523    }
524    out.push('\n');
525    for (i, p) in cloud.points.iter().enumerate() {
526        let _ = write!(out, "{} {} {}", p.x, p.y, p.z);
527        if let Some(normals) = &cloud.normals {
528            let n = normals[i];
529            let _ = write!(out, " {} {} {}", n.x, n.y, n.z);
530        }
531        if let Some(colours) = &cloud.colours {
532            let [r, g, b] = [0, 1, 2].map(|k| channel_byte(colours[i][k]));
533            let _ = write!(out, " {r} {g} {b}");
534        }
535        for s in &cloud.scalars {
536            let _ = write!(out, " {}", s.values[i]);
537        }
538        out.push('\n');
539    }
540    Ok(out)
541}
542
543/// A colour channel in `[0, 1]` as the nearest of 0 to 255.
544pub(crate) fn channel_byte(value: f64) -> u8 {
545    #[allow(
546        clippy::cast_possible_truncation,
547        clippy::cast_sign_loss,
548        reason = "clamped into 0..=255 first"
549    )]
550    let byte = (value.clamp(0.0, 1.0) * 255.0).round() as u8;
551    byte
552}
553
554#[cfg(test)]
555mod tests {
556    #![allow(clippy::unwrap_used, reason = "test code")]
557    use super::*;
558
559    fn sample_cloud() -> PointCloud {
560        PointCloud {
561            points: vec![Point::new(0.1, -2.5, 1e-7), Point::new(3.0, 4.0, 5.0)],
562            normals: Some(vec![Vector::Z, Vector::new(0.6, 0.0, -0.8)]),
563            colours: Some(vec![[1.0, 0.0, 51.0 / 255.0, 1.0], [0.0, 1.0, 0.0, 1.0]]),
564            scalars: vec![Scalar {
565                name: "deviation".into(),
566                values: vec![0.012_345_678_9, -0.25],
567            }],
568        }
569    }
570
571    #[test]
572    fn an_xyz_cloud_round_trips_with_its_header() {
573        let cloud = sample_cloud();
574        let text = write_xyz(&cloud).unwrap();
575        assert!(text.starts_with("# x y z nx ny nz red green blue deviation\n"));
576        assert_eq!(read_xyz(&text, None).unwrap(), cloud);
577    }
578
579    #[test]
580    fn unnamed_columns_are_read_by_their_count() {
581        let rgb = read_xyz("1,2,3,255,0,10\n4,5,6,0,128,0\n", None).unwrap();
582        assert_eq!(rgb.colours.unwrap()[0], [1.0, 0.0, 10.0 / 255.0, 1.0]);
583        assert!(rgb.normals.is_none());
584        let normal = read_xyz("1 2 3 0 0 1\n4 5 6 0 1 0\n", None).unwrap();
585        assert_eq!(normal.normals.unwrap()[1], Vector::Y);
586        let intensity = read_xyz("1 2 3 -40\n", None).unwrap();
587        assert_eq!(intensity.scalar("intensity"), Some(&[-40.0][..]));
588        let nine = read_xyz("0 0 0 0.6 0 0.8 10 20 30\n", None).unwrap();
589        assert_eq!(nine.normals.unwrap()[0], Vector::new(0.6, 0.0, 0.8));
590        assert_eq!(nine.colours.unwrap()[0][1], 20.0 / 255.0);
591        assert!(read_xyz("1 2 3 4 5\n", None).is_err());
592        assert!(read_xyz("0 0 0 10 20 30 10 20 30\n", None).is_err());
593        let named = read_xyz(
594            "1 2 3 9 7\n",
595            Some(&[
596                XyzColumn::Skip,
597                XyzColumn::X,
598                XyzColumn::Y,
599                XyzColumn::Z,
600                XyzColumn::Scalar("quality".into()),
601            ]),
602        )
603        .unwrap();
604        assert_eq!(named.points[0], Point::new(2.0, 3.0, 9.0));
605        assert_eq!(named.scalar("quality"), Some(&[7.0][..]));
606    }
607
608    #[test]
609    fn a_pts_reads_every_scan_and_holds_each_to_its_count() {
610        let text = "2\n0 0 0 -100 255 0 0\n1 0 0 200 0 255 0\n1\n0 1 0 0 0 0 255\n";
611        let cloud = read_pts(text, None).unwrap();
612        assert_eq!(cloud.len(), 3);
613        assert_eq!(cloud.scalar("intensity"), Some(&[-100.0, 200.0, 0.0][..]));
614        assert_eq!(cloud.colours.unwrap()[2], [0.0, 0.0, 1.0, 1.0]);
615        assert!(read_pts("3\n0 0 0\n1 1 1\n", None).is_err());
616        assert!(read_pts("1\n0 0 0\n1 1 1\n", None).is_err());
617        assert!(read_pts("0 0 0\n", None).is_err());
618    }
619
620    #[test]
621    fn garbage_and_truncation_are_refused_by_name() {
622        for bad in [
623            "",
624            "# x y z\n",
625            "1 2 nan\n",
626            "1 2 3\n4 5\n",
627            "1 2 3\n4 5 6 7\n",
628            "1 2 x\n",
629            "1 2 3 0 300 0 1\n",
630            "\u{0}\u{1}\u{2}",
631        ] {
632            assert!(read_xyz(bad, None).is_err(), "{bad:?}");
633        }
634        let text = write_xyz(&sample_cloud()).unwrap();
635        for cut in 1..text.len() {
636            // Every prefix either reads or is refused; none panics.
637            let _ = read_xyz(&text[..cut], None);
638        }
639        let mut bad_name = sample_cloud();
640        bad_name.scalars[0].name = "red".into();
641        assert!(write_xyz(&bad_name).is_err());
642        let mut short = sample_cloud();
643        short.normals = Some(vec![Vector::Z]);
644        assert!(write_xyz(&short).is_err());
645    }
646}
ogeom_io/cloud.rs

ogeom_io/
cloud.rs
