copc_writer/

writer.rs

use std::fs::File;
use std::io::{BufReader, BufWriter, Cursor, Seek, SeekFrom, Write};
use std::path::Path;

use byteorder::{LittleEndian, ReadBytesExt, WriteBytesExt};
use copc_core::{
    Bounds, CancelCheck, ColumnData, CopcInfo, Entry, Error, LasColumnBatch, LasDimension,
    LasPointRecord, NeverCancel, Result, StreamingLayout, VoxelKey, HIERARCHY_ENTRY_BYTES,
};
use las::{point::Format as LasFormat, raw, Color, Read as _};
use laz::{LasZipCompressor, LazVlrBuilder};
use tempfile::{NamedTempFile, TempPath};

use crate::spill::{SpillReader, SpillWriter};

const CANCEL_POLL_STRIDE: usize = 4_096;
const HIERARCHY_PAGE_MAX_ENTRIES: usize = 4_096;
const INDEX_RECORD_BYTES: u64 = 4;
const LAS_14_SCAN_ANGLE_SCALE: f32 = 0.006;
const LASZIP_VLR_USER_ID: &str = "laszip encoded";
const LASZIP_VLR_RECORD_ID: u16 = 22204;

/// Normalized point fields consumed by the COPC writer.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CopcPointFields {
    pub x: f64,
    pub y: f64,
    pub z: f64,
    pub intensity: u16,
    pub return_number: u8,
    pub number_of_returns: u8,
    pub synthetic: u8,
    pub key_point: u8,
    pub withheld: u8,
    pub overlap: u8,
    pub scan_channel: u8,
    pub scan_direction_flag: u8,
    pub edge_of_flight_line: u8,
    pub classification: u8,
    pub user_data: u8,
    /// Scan angle in degrees; encoded as LAS 1.4 scaled scan angle on write.
    pub scan_angle: f32,
    pub point_source_id: u16,
    pub gps_time: f64,
    pub red: u16,
    pub green: u16,
    pub blue: u16,
}

/// Abstract point-data source for COPC emission.
pub trait CopcPointSource {
    fn len(&self) -> usize;
    fn xyz(&self, index: usize) -> (f64, f64, f64);
    fn fields(&self, index: usize) -> Result<CopcPointFields>;

    fn is_empty(&self) -> bool {
        self.len() == 0
    }
}

/// COPC writer source backed directly by a neutral LAS column batch.
pub struct ColumnBatchSource<'a> {
    batch: &'a LasColumnBatch,
    x: &'a [f64],
    y: &'a [f64],
    z: &'a [f64],
    intensity: Option<&'a [u16]>,
    return_number: Option<&'a [u8]>,
    number_of_returns: Option<&'a [u8]>,
    synthetic: Option<&'a [bool]>,
    key_point: Option<&'a [bool]>,
    withheld: Option<&'a [bool]>,
    overlap: Option<&'a [bool]>,
    scan_channel: Option<&'a [u8]>,
    scan_direction_flag: Option<&'a [bool]>,
    edge_of_flight_line: Option<&'a [bool]>,
    classification: Option<&'a [u8]>,
    user_data: Option<&'a [u8]>,
    scan_angle_rank: Option<&'a [i16]>,
    point_source_id: Option<&'a [u16]>,
    gps_time: Option<&'a [f64]>,
    red: Option<&'a [u16]>,
    green: Option<&'a [u16]>,
    blue: Option<&'a [u16]>,
}

impl<'a> ColumnBatchSource<'a> {
    pub fn new(batch: &'a LasColumnBatch) -> Result<Self> {
        batch.validate()?;
        validate_column_batch_writer_support(batch)?;

        let x = required_f64_column(batch, LasDimension::X)?;
        let y = required_f64_column(batch, LasDimension::Y)?;
        let z = required_f64_column(batch, LasDimension::Z)?;
        let red = optional_u16_column(batch, LasDimension::Red)?;
        let green = optional_u16_column(batch, LasDimension::Green)?;
        let blue = optional_u16_column(batch, LasDimension::Blue)?;
        validate_color_columns(red, green, blue)?;

        Ok(Self {
            batch,
            x,
            y,
            z,
            intensity: optional_u16_column(batch, LasDimension::Intensity)?,
            return_number: optional_u8_column(batch, LasDimension::ReturnNumber)?,
            number_of_returns: optional_u8_column(batch, LasDimension::NumberOfReturns)?,
            synthetic: optional_bool_column(batch, LasDimension::Synthetic)?,
            key_point: optional_bool_column(batch, LasDimension::KeyPoint)?,
            withheld: optional_bool_column(batch, LasDimension::Withheld)?,
            overlap: optional_bool_column(batch, LasDimension::Overlap)?,
            scan_channel: optional_u8_column(batch, LasDimension::ScanChannel)?,
            scan_direction_flag: optional_bool_column(batch, LasDimension::ScanDirectionFlag)?,
            edge_of_flight_line: optional_bool_column(batch, LasDimension::EdgeOfFlightLine)?,
            classification: optional_u8_column(batch, LasDimension::Classification)?,
            user_data: optional_u8_column(batch, LasDimension::UserData)?,
            scan_angle_rank: optional_i16_column(batch, LasDimension::ScanAngleRank)?,
            point_source_id: optional_u16_column(batch, LasDimension::PointSourceId)?,
            gps_time: optional_f64_column(batch, LasDimension::GpsTime)?,
            red,
            green,
            blue,
        })
    }

    pub fn batch(&self) -> &LasColumnBatch {
        self.batch
    }

    pub fn has_color(&self) -> bool {
        self.red.is_some() && self.green.is_some() && self.blue.is_some()
    }

    pub fn bounds(&self) -> Result<Bounds> {
        if self.is_empty() {
            return Err(Error::InvalidInput(
                "cannot compute bounds for empty column batch".into(),
            ));
        }
        let mut bounds = Bounds::point(self.x[0], self.y[0], self.z[0]);
        for index in 1..self.len() {
            bounds.extend(self.x[index], self.y[index], self.z[index]);
        }
        Ok(bounds)
    }
}

impl CopcPointSource for ColumnBatchSource<'_> {
    fn len(&self) -> usize {
        self.batch.len()
    }

    #[inline]
    fn xyz(&self, index: usize) -> (f64, f64, f64) {
        (self.x[index], self.y[index], self.z[index])
    }

    fn fields(&self, index: usize) -> Result<CopcPointFields> {
        Ok(CopcPointFields {
            x: self.x[index],
            y: self.y[index],
            z: self.z[index],
            intensity: at_u16(self.intensity, index),
            return_number: at_u8(self.return_number, index),
            number_of_returns: at_u8(self.number_of_returns, index),
            synthetic: at_bool_u8(self.synthetic, index),
            key_point: at_bool_u8(self.key_point, index),
            withheld: at_bool_u8(self.withheld, index),
            overlap: at_bool_u8(self.overlap, index),
            scan_channel: at_u8(self.scan_channel, index),
            scan_direction_flag: at_bool_u8(self.scan_direction_flag, index),
            edge_of_flight_line: at_bool_u8(self.edge_of_flight_line, index),
            classification: at_u8(self.classification, index),
            user_data: at_u8(self.user_data, index),
            scan_angle: self
                .scan_angle_rank
                .map(|column| column[index] as f32 * 90.0 / 180.0)
                .unwrap_or(0.0),
            point_source_id: at_u16(self.point_source_id, index),
            gps_time: self.gps_time.map(|column| column[index]).unwrap_or(0.0),
            red: at_u16(self.red, index),
            green: at_u16(self.green, index),
            blue: at_u16(self.blue, index),
        })
    }
}

fn at_bool_u8(column: Option<&[bool]>, index: usize) -> u8 {
    column.map(|values| u8::from(values[index])).unwrap_or(0)
}

fn at_u8(column: Option<&[u8]>, index: usize) -> u8 {
    column.map(|values| values[index]).unwrap_or(0)
}

fn at_u16(column: Option<&[u16]>, index: usize) -> u16 {
    column.map(|values| values[index]).unwrap_or(0)
}

fn validate_column_batch_writer_support(batch: &LasColumnBatch) -> Result<()> {
    let unsupported: Vec<_> = batch
        .columns
        .iter()
        .filter_map(|(spec, _)| match spec.dimension {
            LasDimension::Nir => Some("NIR point data"),
            LasDimension::WaveformPacketDescriptorIndex
            | LasDimension::WaveformPacketByteOffset
            | LasDimension::WaveformPacketSize
            | LasDimension::WavePacketReturnPointWaveformLocation => Some("waveform point data"),
            LasDimension::ExtraBytes => Some("extra point bytes"),
            _ => None,
        })
        .collect();
    if unsupported.is_empty() {
        Ok(())
    } else {
        Err(Error::Unsupported(format!(
            "COPC writer cannot preserve {}",
            unsupported.join(", ")
        )))
    }
}

fn validate_color_columns(
    red: Option<&[u16]>,
    green: Option<&[u16]>,
    blue: Option<&[u16]>,
) -> Result<()> {
    let present =
        usize::from(red.is_some()) + usize::from(green.is_some()) + usize::from(blue.is_some());
    if present == 0 || present == 3 {
        Ok(())
    } else {
        Err(Error::InvalidInput(
            "Red, Green, and Blue columns must be supplied together".into(),
        ))
    }
}

fn required_f64_column(batch: &LasColumnBatch, dimension: LasDimension) -> Result<&[f64]> {
    match batch.column(dimension) {
        Some(ColumnData::F64(values)) => Ok(values),
        Some(other) => Err(unexpected_column_type(dimension, "F64", other)),
        None => Err(Error::InvalidInput(format!(
            "ColumnBatchSource requires {dimension:?} column"
        ))),
    }
}

fn optional_f64_column(batch: &LasColumnBatch, dimension: LasDimension) -> Result<Option<&[f64]>> {
    match batch.column(dimension) {
        Some(ColumnData::F64(values)) => Ok(Some(values)),
        Some(other) => Err(unexpected_column_type(dimension, "F64", other)),
        None => Ok(None),
    }
}

fn optional_i16_column(batch: &LasColumnBatch, dimension: LasDimension) -> Result<Option<&[i16]>> {
    match batch.column(dimension) {
        Some(ColumnData::I16(values)) => Ok(Some(values)),
        Some(other) => Err(unexpected_column_type(dimension, "I16", other)),
        None => Ok(None),
    }
}

fn optional_u16_column(batch: &LasColumnBatch, dimension: LasDimension) -> Result<Option<&[u16]>> {
    match batch.column(dimension) {
        Some(ColumnData::U16(values)) => Ok(Some(values)),
        Some(other) => Err(unexpected_column_type(dimension, "U16", other)),
        None => Ok(None),
    }
}

fn optional_u8_column(batch: &LasColumnBatch, dimension: LasDimension) -> Result<Option<&[u8]>> {
    match batch.column(dimension) {
        Some(ColumnData::U8(values)) => Ok(Some(values)),
        Some(other) => Err(unexpected_column_type(dimension, "U8", other)),
        None => Ok(None),
    }
}

fn optional_bool_column(
    batch: &LasColumnBatch,
    dimension: LasDimension,
) -> Result<Option<&[bool]>> {
    match batch.column(dimension) {
        Some(ColumnData::Bool(values)) => Ok(Some(values)),
        Some(other) => Err(unexpected_column_type(dimension, "Bool", other)),
        None => Ok(None),
    }
}

fn unexpected_column_type(dimension: LasDimension, expected: &str, actual: &ColumnData) -> Error {
    Error::InvalidInput(format!(
        "{dimension:?} column must be {expected}, found {:?}",
        actual.scalar()
    ))
}

struct SpillSource<'a> {
    reader: &'a SpillReader,
}

impl CopcPointSource for SpillSource<'_> {
    fn len(&self) -> usize {
        self.reader.len()
    }

    #[inline]
    fn xyz(&self, index: usize) -> (f64, f64, f64) {
        self.reader.xyz_at(index)
    }

    fn fields(&self, index: usize) -> Result<CopcPointFields> {
        let record = self.reader.record_at(index)?;
        Ok(CopcPointFields {
            x: record.x,
            y: record.y,
            z: record.z,
            intensity: record.intensity,
            return_number: record.return_number,
            number_of_returns: record.number_of_returns,
            synthetic: u8::from(record.synthetic),
            key_point: u8::from(record.key_point),
            withheld: u8::from(record.withheld),
            overlap: u8::from(record.overlap),
            scan_channel: record.scan_channel,
            scan_direction_flag: u8::from(record.scan_direction_flag),
            edge_of_flight_line: u8::from(record.edge_of_flight_line),
            classification: record.classification,
            user_data: record.user_data,
            scan_angle: record.scan_angle,
            point_source_id: record.point_source_id,
            gps_time: record.gps_time,
            red: record.red,
            green: record.green,
            blue: record.blue,
        })
    }
}

#[derive(Clone, Debug)]
struct OutputLasMetadata {
    file_source_id: u16,
    global_encoding: u16,
    guid: [u8; 16],
    system_identifier: String,
    generating_software: String,
    creation_day_of_year: u16,
    creation_year: u16,
    scale: (f64, f64, f64),
    offset: Option<(f64, f64, f64)>,
}

impl Default for OutputLasMetadata {
    fn default() -> Self {
        Self {
            file_source_id: 0,
            global_encoding: 0,
            guid: [0; 16],
            system_identifier: "copc-rust".to_string(),
            generating_software: "copc-writer".to_string(),
            creation_day_of_year: 0,
            creation_year: 2026,
            scale: (0.001, 0.001, 0.001),
            offset: None,
        }
    }
}

impl OutputLasMetadata {
    fn from_las_header(header: &las::Header) -> Self {
        let mut global_encoding = u16::from(header.gps_time_type());
        if header.has_synthetic_return_numbers() {
            global_encoding |= 8;
        }
        let transforms = header.transforms();
        let (creation_day_of_year, creation_year) = header
            .date()
            .map(|date| {
                let year = date.format("%Y").to_string().parse().unwrap_or(0);
                let day = date.format("%j").to_string().parse().unwrap_or(0);
                (day, year)
            })
            .unwrap_or((0, 0));

        Self {
            file_source_id: header.file_source_id(),
            global_encoding,
            guid: *header.guid().as_bytes(),
            system_identifier: header.system_identifier().to_string(),
            generating_software: header.generating_software().to_string(),
            creation_day_of_year,
            creation_year,
            scale: (transforms.x.scale, transforms.y.scale, transforms.z.scale),
            offset: Some((
                transforms.x.offset,
                transforms.y.offset,
                transforms.z.offset,
            )),
        }
    }
}

#[derive(Clone, Copy, Debug, PartialEq)]
struct PointStats {
    gpstime_min: f64,
    gpstime_max: f64,
    extended_return_counts: [u64; 15],
}

impl PointStats {
    fn new() -> Self {
        Self {
            gpstime_min: f64::INFINITY,
            gpstime_max: f64::NEG_INFINITY,
            extended_return_counts: [0; 15],
        }
    }

    fn record(&mut self, index: usize, fields: &CopcPointFields) -> Result<()> {
        validate_finite_value(&format!("point {index} GPS time"), fields.gps_time)?;
        self.gpstime_min = self.gpstime_min.min(fields.gps_time);
        self.gpstime_max = self.gpstime_max.max(fields.gps_time);
        if (1..=15).contains(&fields.return_number) {
            self.extended_return_counts[usize::from(fields.return_number - 1)] += 1;
        }
        Ok(())
    }
}

#[derive(Debug, Clone, Copy)]
pub struct CopcWriterParams {
    pub max_points_per_node: u32,
    pub max_depth: u32,
}

impl Default for CopcWriterParams {
    fn default() -> Self {
        Self {
            max_points_per_node: 100_000,
            max_depth: 8,
        }
    }
}

pub fn write_source<S: CopcPointSource>(
    path: &Path,
    source: &S,
    has_color: bool,
    bounds: Bounds,
    params: &CopcWriterParams,
) -> Result<()> {
    write_source_with_cancel(path, source, has_color, bounds, params, &NeverCancel)
}

pub fn write_source_with_cancel<S: CopcPointSource>(
    path: &Path,
    source: &S,
    has_color: bool,
    bounds: Bounds,
    params: &CopcWriterParams,
    cancel: &dyn CancelCheck,
) -> Result<()> {
    cancel.check()?;
    if source.is_empty() {
        return Err(Error::InvalidInput(
            "cannot write empty cloud to COPC".into(),
        ));
    }
    write_copc_inner(
        path,
        source,
        has_color,
        bounds,
        params,
        cancel,
        &OutputLasMetadata::default(),
    )
}

pub fn write_streaming_with_cancel<I>(
    path: &Path,
    layout: StreamingLayout,
    points: I,
    params: &CopcWriterParams,
    spill_dir: &Path,
    cancel: &dyn CancelCheck,
) -> Result<()>
where
    I: IntoIterator<Item = Result<LasPointRecord>>,
{
    cancel.check()?;
    validate_streaming_layout_supported(&layout)?;
    let mut spill = SpillWriter::create(spill_dir, layout)?;
    for (index, item) in points.into_iter().enumerate() {
        if index % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        let record = item?;
        validate_record_coordinates(&record, index)?;
        spill.push(&record)?;
    }
    cancel.check()?;
    let reader = spill.finalize()?;
    write_copc_from_spill(path, reader, params, cancel, &OutputLasMetadata::default())
}

pub fn convert_las_to_copc_streaming(
    las_path: &Path,
    copc_path: &Path,
    params: &CopcWriterParams,
    spill_dir: &Path,
    cancel: &dyn CancelCheck,
) -> Result<()> {
    cancel.check()?;
    let mut reader = las::Reader::from_path(las_path).map_err(|e| Error::Las(e.to_string()))?;
    validate_las_conversion_supported(reader.header())?;
    let output_metadata = OutputLasMetadata::from_las_header(reader.header());
    let layout = StreamingLayout::from_las_format(*reader.header().point_format());
    let mut spill = SpillWriter::create(spill_dir, layout)?;
    for (index, result) in reader.points().enumerate() {
        if index % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        let point = result.map_err(|e| Error::Las(e.to_string()))?;
        let record = LasPointRecord::from_las_point(&point);
        validate_record_coordinates(&record, index)?;
        spill.push(&record)?;
    }
    cancel.check()?;
    let reader = spill.finalize()?;
    write_copc_from_spill(copc_path, reader, params, cancel, &output_metadata)
}

fn validate_streaming_layout_supported(layout: &StreamingLayout) -> Result<()> {
    let mut unsupported = Vec::new();
    if layout.has_nir {
        unsupported.push("NIR point data".to_string());
    }
    if layout.has_waveform {
        unsupported.push("waveform point data".to_string());
    }
    if unsupported.is_empty() {
        Ok(())
    } else {
        Err(Error::Unsupported(format!(
            "COPC writer cannot preserve {}",
            unsupported.join(", ")
        )))
    }
}

fn validate_las_conversion_supported(header: &las::Header) -> Result<()> {
    let mut unsupported = Vec::new();
    let format = header.point_format();
    if format.has_nir {
        unsupported.push("NIR point data".to_string());
    }
    if format.has_waveform {
        unsupported.push("waveform point data".to_string());
    }
    if format.extra_bytes > 0 {
        unsupported.push(format!("{} extra point byte(s)", format.extra_bytes));
    }
    let unsupported_vlr_count = header
        .vlrs()
        .iter()
        .filter(|vlr| !is_laszip_vlr(vlr))
        .count();
    if unsupported_vlr_count > 0 {
        unsupported.push(format!("{unsupported_vlr_count} VLR(s)"));
    }
    if !header.evlrs().is_empty() {
        unsupported.push(format!("{} EVLR(s)", header.evlrs().len()));
    }
    if !header.padding().is_empty() {
        unsupported.push(format!("{} header padding byte(s)", header.padding().len()));
    }
    if !header.vlr_padding().is_empty() {
        unsupported.push(format!(
            "{} VLR padding byte(s)",
            header.vlr_padding().len()
        ));
    }
    if !header.point_padding().is_empty() {
        unsupported.push(format!(
            "{} point padding byte(s)",
            header.point_padding().len()
        ));
    }

    if unsupported.is_empty() {
        Ok(())
    } else {
        Err(Error::Unsupported(format!(
            "LAS-to-COPC streaming conversion cannot preserve {}",
            unsupported.join(", ")
        )))
    }
}

fn is_laszip_vlr(vlr: &las::Vlr) -> bool {
    vlr.user_id == LASZIP_VLR_USER_ID && vlr.record_id == LASZIP_VLR_RECORD_ID
}

fn validate_record_coordinates(record: &LasPointRecord, index: usize) -> Result<()> {
    validate_xyz_finite(index, record.x, record.y, record.z)
}

fn validate_coordinate_inputs<S: CopcPointSource>(
    source: &S,
    bounds: Bounds,
    scale: (f64, f64, f64),
    offset: (f64, f64, f64),
    cancel: &dyn CancelCheck,
) -> Result<PointStats> {
    validate_bounds(bounds)?;
    validate_transform(scale, offset)?;
    let mut stats = PointStats::new();
    for index in 0..source.len() {
        if index % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        let (x, y, z) = source.xyz(index);
        validate_xyz_finite(index, x, y, z)?;
        quantize_xyz(index, x, y, z, scale, offset)?;

        let fields = source.fields(index)?;
        validate_xyz_finite(index, fields.x, fields.y, fields.z)?;
        quantize_xyz(index, fields.x, fields.y, fields.z, scale, offset)?;
        validate_scan_angle(index, fields.scan_angle)?;
        validate_point_flags(index, &fields)?;
        stats.record(index, &fields)?;
    }
    Ok(stats)
}

fn validate_point_flags(index: usize, fields: &CopcPointFields) -> Result<()> {
    validate_point_field_range(index, "return_number", fields.return_number, 0, 15)?;
    validate_point_field_range(index, "number_of_returns", fields.number_of_returns, 0, 15)?;
    validate_point_field_range(index, "synthetic", fields.synthetic, 0, 1)?;
    validate_point_field_range(index, "key_point", fields.key_point, 0, 1)?;
    validate_point_field_range(index, "withheld", fields.withheld, 0, 1)?;
    validate_point_field_range(index, "overlap", fields.overlap, 0, 1)?;
    validate_point_field_range(index, "scan_channel", fields.scan_channel, 0, 3)?;
    validate_point_field_range(
        index,
        "scan_direction_flag",
        fields.scan_direction_flag,
        0,
        1,
    )?;
    validate_point_field_range(
        index,
        "edge_of_flight_line",
        fields.edge_of_flight_line,
        0,
        1,
    )
}

fn validate_scan_angle(index: usize, value: f32) -> Result<()> {
    if !value.is_finite() {
        return Err(Error::InvalidInput(format!(
            "point {index} scan_angle must be finite, got {value}"
        )));
    }
    let scaled = value / LAS_14_SCAN_ANGLE_SCALE;
    if scaled < f32::from(i16::MIN) || scaled > f32::from(i16::MAX) {
        return Err(Error::InvalidInput(format!(
            "point {index} scan_angle {value} encodes to {scaled}, outside LAS i16 range"
        )));
    }
    Ok(())
}

fn validate_point_field_range(index: usize, name: &str, value: u8, min: u8, max: u8) -> Result<()> {
    if (min..=max).contains(&value) {
        Ok(())
    } else {
        Err(Error::InvalidInput(format!(
            "point {index} {name} must be in {min}..={max}, got {value}"
        )))
    }
}

fn validate_bounds(bounds: Bounds) -> Result<()> {
    validate_finite_value("bounds min x", bounds.min.0)?;
    validate_finite_value("bounds min y", bounds.min.1)?;
    validate_finite_value("bounds min z", bounds.min.2)?;
    validate_finite_value("bounds max x", bounds.max.0)?;
    validate_finite_value("bounds max y", bounds.max.1)?;
    validate_finite_value("bounds max z", bounds.max.2)?;
    for (axis, min, max) in [
        ("x", bounds.min.0, bounds.max.0),
        ("y", bounds.min.1, bounds.max.1),
        ("z", bounds.min.2, bounds.max.2),
    ] {
        if min > max {
            return Err(Error::InvalidInput(format!(
                "bounds {axis} min {min} exceeds max {max}"
            )));
        }
        validate_finite_value(&format!("bounds {axis} span"), max - min)?;
    }
    Ok(())
}

fn validate_transform(scale: (f64, f64, f64), offset: (f64, f64, f64)) -> Result<()> {
    for (axis, value) in [("x", scale.0), ("y", scale.1), ("z", scale.2)] {
        if !value.is_finite() || value <= 0.0 {
            return Err(Error::InvalidInput(format!(
                "LAS {axis} scale must be finite and positive, got {value}"
            )));
        }
    }
    validate_finite_value("LAS x offset", offset.0)?;
    validate_finite_value("LAS y offset", offset.1)?;
    validate_finite_value("LAS z offset", offset.2)?;
    Ok(())
}

fn validate_xyz_finite(index: usize, x: f64, y: f64, z: f64) -> Result<()> {
    validate_point_axis_finite(index, "x", x)?;
    validate_point_axis_finite(index, "y", y)?;
    validate_point_axis_finite(index, "z", z)
}

fn validate_point_axis_finite(index: usize, axis: &str, value: f64) -> Result<()> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(Error::InvalidInput(format!(
            "point {index} {axis} coordinate must be finite, got {value}"
        )))
    }
}

fn validate_finite_value(name: &str, value: f64) -> Result<()> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(Error::InvalidInput(format!(
            "{name} must be finite, got {value}"
        )))
    }
}

fn quantize_xyz(
    index: usize,
    x: f64,
    y: f64,
    z: f64,
    scale: (f64, f64, f64),
    offset: (f64, f64, f64),
) -> Result<(i32, i32, i32)> {
    Ok((
        quantize_axis(index, "x", x, scale.0, offset.0)?,
        quantize_axis(index, "y", y, scale.1, offset.1)?,
        quantize_axis(index, "z", z, scale.2, offset.2)?,
    ))
}

fn quantize_axis(index: usize, axis: &str, value: f64, scale: f64, offset: f64) -> Result<i32> {
    let scaled = ((value - offset) / scale).round();
    if !scaled.is_finite() {
        return Err(Error::InvalidInput(format!(
            "point {index} {axis} coordinate cannot be encoded with scale {scale} and offset {offset}"
        )));
    }
    if scaled < f64::from(i32::MIN) || scaled > f64::from(i32::MAX) {
        return Err(Error::InvalidInput(format!(
            "point {index} {axis} coordinate {value} encodes to {scaled}, outside LAS i32 range"
        )));
    }
    Ok(scaled as i32)
}

fn write_copc_from_spill(
    path: &Path,
    reader: SpillReader,
    params: &CopcWriterParams,
    cancel: &dyn CancelCheck,
    metadata: &OutputLasMetadata,
) -> Result<()> {
    cancel.check()?;
    validate_streaming_layout_supported(&reader.layout())?;
    if reader.is_empty() {
        return Err(Error::InvalidInput(
            "cannot write empty cloud to COPC".into(),
        ));
    }
    let has_color = reader.layout().has_color;
    let bounds = reader.bounds();
    let source = SpillSource { reader: &reader };
    write_copc_inner(path, &source, has_color, bounds, params, cancel, metadata)
}

fn write_copc_inner<S: CopcPointSource>(
    path: &Path,
    source: &S,
    has_color: bool,
    bounds: Bounds,
    params: &CopcWriterParams,
    cancel: &dyn CancelCheck,
    metadata: &OutputLasMetadata,
) -> Result<()> {
    cancel.check()?;
    let point_format_id = if has_color { 7u8 } else { 6u8 };
    let point_format =
        LasFormat::new(point_format_id).map_err(|e| Error::Las(format!("point format: {e}")))?;
    let point_record_length = point_format.len();

    let (scale_x, scale_y, scale_z) = metadata.scale;
    let (offset_x, offset_y, offset_z) =
        metadata
            .offset
            .unwrap_or((bounds.min.0, bounds.min.1, bounds.min.2));
    let point_stats = validate_coordinate_inputs(
        source,
        bounds,
        (scale_x, scale_y, scale_z),
        (offset_x, offset_y, offset_z),
        cancel,
    )?;
    let (center, halfsize) = cube_from_bounds(&bounds);

    let lod_index = build_lod_index(source, center, halfsize, params, cancel)?;
    cancel.check()?;

    let var_vlr = LazVlrBuilder::default()
        .with_point_format(point_format_id, 0)
        .map_err(|e| Error::Las(format!("laz items: {e}")))?
        .with_variable_chunk_size()
        .build();
    let mut var_vlr_bytes = Vec::new();
    var_vlr
        .write_to(&mut var_vlr_bytes)
        .map_err(|e| Error::Las(format!("variable chunk LAZ VLR: {e}")))?;

    let copc_info_vlr_size = 160u16;
    let las_header_size = 375u32;
    let total_vlr_bytes =
        (54u32 + u32::from(copc_info_vlr_size)) + (54u32 + var_vlr_bytes.len() as u32);
    let offset_to_point_data = las_header_size + total_vlr_bytes;

    let file = File::create(path).map_err(|e| Error::io("create COPC file", e))?;
    let mut writer = BufWriter::new(file);

    let header = LasHeader {
        point_data_format: point_format_id | 0x80,
        point_record_length,
        offset_to_point_data,
        number_of_vlrs: 2,
        file_source_id: metadata.file_source_id,
        global_encoding: metadata.global_encoding,
        guid: metadata.guid,
        system_identifier: metadata.system_identifier.clone(),
        generating_software: metadata.generating_software.clone(),
        creation_day_of_year: metadata.creation_day_of_year,
        creation_year: metadata.creation_year,
        scale: (scale_x, scale_y, scale_z),
        offset: (offset_x, offset_y, offset_z),
        bounds,
        legacy_point_count: 0,
        total_point_count: source.len() as u64,
        offset_to_first_evlr: 0,
        number_of_evlrs: 1,
        extended_return_counts: point_stats.extended_return_counts,
    };
    header.write(&mut writer)?;

    write_vlr_header(&mut writer, "copc", 1, copc_info_vlr_size, "COPC info")?;
    let copc_info_payload_start = writer
        .stream_position()
        .map_err(|e| Error::io("record COPC info payload offset", e))?;
    writer
        .write_all(&[0u8; 160])
        .map_err(|e| Error::io("write COPC info placeholder", e))?;

    write_vlr_header(
        &mut writer,
        LASZIP_VLR_USER_ID,
        LASZIP_VLR_RECORD_ID,
        var_vlr_bytes.len() as u16,
        "http://laszip.org",
    )?;
    writer
        .write_all(&var_vlr_bytes)
        .map_err(|e| Error::io("write LAZ VLR", e))?;

    let point_data_actual_start = writer
        .stream_position()
        .map_err(|e| Error::io("record point data offset", e))?;
    if point_data_actual_start as u32 != offset_to_point_data {
        return Err(Error::InvalidInput(format!(
            "VLR size accounting mismatch: at {point_data_actual_start}, expected {offset_to_point_data}"
        )));
    }

    let mut compressor = LasZipCompressor::new(&mut writer, var_vlr.clone())
        .map_err(|e| Error::Las(format!("compressor: {e}")))?;
    let mut hierarchy: Vec<Entry> = Vec::with_capacity(lod_index.nodes.len());
    let order_path: &Path = lod_index.order_path.as_ref();
    let mut index_reader =
        BufReader::new(File::open(order_path).map_err(|e| Error::io("open LOD order", e))?);
    let mut point_buf = vec![0u8; point_record_length as usize];
    let mut chunk_start_file_offset = compressor
        .get_mut()
        .stream_position()
        .map_err(|e| Error::io("record chunk start", e))?;
    chunk_start_file_offset += 8;

    for node in &lod_index.nodes {
        cancel.check()?;
        index_reader
            .seek(SeekFrom::Start(node.start))
            .map_err(|e| Error::io("seek LOD order", e))?;
        for point_index in 0..node.count {
            if point_index % CANCEL_POLL_STRIDE == 0 {
                cancel.check()?;
            }
            let source_index = index_reader
                .read_u32::<LittleEndian>()
                .map_err(|e| Error::io("read LOD order", e))?
                as usize;
            let fields = source.fields(source_index as usize)?;
            encode_point_record(
                &mut point_buf,
                &fields,
                (scale_x, scale_y, scale_z),
                (offset_x, offset_y, offset_z),
                source_index as usize,
                &point_format,
            )?;
            compressor
                .compress_one(&point_buf)
                .map_err(|e| Error::Las(format!("compress point: {e}")))?;
        }
        compressor
            .finish_current_chunk()
            .map_err(|e| Error::Las(format!("finish chunk: {e}")))?;
        let after = compressor
            .get_mut()
            .stream_position()
            .map_err(|e| Error::io("record chunk end", e))?;
        let byte_size = i32::try_from(after - chunk_start_file_offset)
            .map_err(|_| Error::InvalidInput("LAZ chunk exceeds COPC i32 byte size".into()))?;
        let point_count = i32::try_from(node.count)
            .map_err(|_| Error::InvalidInput("node point count exceeds COPC i32 range".into()))?;
        hierarchy.push(Entry {
            key: node.key,
            offset: chunk_start_file_offset,
            byte_size,
            point_count,
        });
        chunk_start_file_offset = after;
    }

    cancel.check()?;
    compressor
        .done()
        .map_err(|e| Error::Las(format!("finish compressor: {e}")))?;
    drop(compressor);

    let evlr_start = writer
        .stream_position()
        .map_err(|e| Error::io("record EVLR start", e))?;
    let root_hier_offset = evlr_start
        .checked_add(60)
        .ok_or_else(|| Error::InvalidInput("hierarchy EVLR offset overflow".into()))?;
    let mut hierarchy_pages = plan_hierarchy_pages(&hierarchy, VoxelKey::root())?;
    let hierarchy_end = assign_hierarchy_page_offsets(&mut hierarchy_pages, root_hier_offset)?;
    let hierarchy_body_size = hierarchy_end
        .checked_sub(root_hier_offset)
        .ok_or_else(|| Error::InvalidInput("hierarchy size overflow".into()))?;
    write_evlr_header(
        &mut writer,
        "copc",
        1000,
        hierarchy_body_size,
        "COPC hierarchy",
    )?;
    let actual_root_hier_offset = writer
        .stream_position()
        .map_err(|e| Error::io("record root hierarchy offset", e))?;
    if actual_root_hier_offset != root_hier_offset {
        return Err(Error::InvalidInput(format!(
            "hierarchy offset accounting mismatch: at {actual_root_hier_offset}, expected {root_hier_offset}"
        )));
    }
    write_hierarchy_page_tree(&mut writer, &hierarchy_pages)?;

    writer
        .seek(SeekFrom::Start(copc_info_payload_start))
        .map_err(|e| Error::io("seek COPC info payload", e))?;
    let info = CopcInfo {
        center,
        halfsize,
        spacing: halfsize / 128.0,
        root_hier_offset,
        root_hier_size: hierarchy_pages.byte_size,
        gpstime_min: point_stats.gpstime_min,
        gpstime_max: point_stats.gpstime_max,
    };
    writer
        .write_all(&info.write_le_bytes())
        .map_err(|e| Error::io("patch COPC info", e))?;

    writer
        .seek(SeekFrom::Start(235))
        .map_err(|e| Error::io("seek first EVLR offset", e))?;
    writer
        .write_u64::<LittleEndian>(evlr_start)
        .map_err(|e| Error::io("patch first EVLR offset", e))?;

    writer
        .flush()
        .map_err(|e| Error::io("flush COPC file", e))?;
    Ok(())
}

#[derive(Debug)]
struct HierarchyPagePlan {
    key: VoxelKey,
    items: Vec<HierarchyPageItem>,
    offset: u64,
    byte_size: u64,
}

#[derive(Debug)]
enum HierarchyPageItem {
    Point(Entry),
    Child(Box<HierarchyPagePlan>),
}

fn plan_hierarchy_pages(entries: &[Entry], key: VoxelKey) -> Result<HierarchyPagePlan> {
    if entries.is_empty() {
        return Err(Error::InvalidInput(
            "cannot write empty hierarchy page".into(),
        ));
    }
    if entries.len() <= HIERARCHY_PAGE_MAX_ENTRIES {
        return Ok(HierarchyPagePlan {
            key,
            items: entries
                .iter()
                .copied()
                .map(HierarchyPageItem::Point)
                .collect(),
            offset: 0,
            byte_size: 0,
        });
    }

    let mut point_entry = None;
    let mut child_entries: [Vec<Entry>; 8] = std::array::from_fn(|_| Vec::new());
    for entry in entries.iter().copied() {
        if entry.key == key {
            point_entry = Some(entry);
            continue;
        }
        let mut matched = false;
        for (octant, child_entries) in child_entries.iter_mut().enumerate() {
            let child_key = key.child(octant as u8);
            if key_contains(child_key, entry.key) {
                child_entries.push(entry);
                matched = true;
                break;
            }
        }
        if !matched {
            return Err(Error::InvalidInput(format!(
                "hierarchy entry {:?} is not under page key {:?}",
                entry.key, key
            )));
        }
    }

    let mut items = Vec::new();
    if let Some(entry) = point_entry {
        items.push(HierarchyPageItem::Point(entry));
    }
    for (octant, child_entries) in child_entries.into_iter().enumerate() {
        if child_entries.is_empty() {
            continue;
        }
        items.push(HierarchyPageItem::Child(Box::new(plan_hierarchy_pages(
            &child_entries,
            key.child(octant as u8),
        )?)));
    }
    if items.len() > HIERARCHY_PAGE_MAX_ENTRIES {
        return Err(Error::InvalidInput(format!(
            "hierarchy page for {:?} has {} entries, max is {}",
            key,
            items.len(),
            HIERARCHY_PAGE_MAX_ENTRIES
        )));
    }
    Ok(HierarchyPagePlan {
        key,
        items,
        offset: 0,
        byte_size: 0,
    })
}

fn assign_hierarchy_page_offsets(page: &mut HierarchyPagePlan, offset: u64) -> Result<u64> {
    page.offset = offset;
    page.byte_size = hierarchy_page_byte_size(page.items.len())?;
    let mut next = offset
        .checked_add(page.byte_size)
        .ok_or_else(|| Error::InvalidInput("hierarchy page offset overflow".into()))?;
    for item in &mut page.items {
        if let HierarchyPageItem::Child(child) = item {
            next = assign_hierarchy_page_offsets(child, next)?;
        }
    }
    Ok(next)
}

fn hierarchy_page_byte_size(entry_count: usize) -> Result<u64> {
    let bytes = entry_count
        .checked_mul(HIERARCHY_ENTRY_BYTES)
        .ok_or_else(|| Error::InvalidInput("hierarchy page size overflow".into()))?;
    u64::try_from(bytes).map_err(|_| Error::InvalidInput("hierarchy page is too large".into()))
}

fn write_hierarchy_page_tree<W: Write + Seek>(
    writer: &mut W,
    page: &HierarchyPagePlan,
) -> Result<()> {
    let position = writer
        .stream_position()
        .map_err(|e| Error::io("record hierarchy page offset", e))?;
    if position != page.offset {
        return Err(Error::InvalidInput(format!(
            "hierarchy page offset mismatch: at {position}, expected {}",
            page.offset
        )));
    }
    let mut entry_buf = [0u8; HIERARCHY_ENTRY_BYTES];
    for item in &page.items {
        hierarchy_page_item_entry(item)?.write_le(&mut entry_buf)?;
        writer
            .write_all(&entry_buf)
            .map_err(|e| Error::io("write hierarchy entry", e))?;
    }
    for item in &page.items {
        if let HierarchyPageItem::Child(child) = item {
            write_hierarchy_page_tree(writer, child)?;
        }
    }
    Ok(())
}

fn hierarchy_page_item_entry(item: &HierarchyPageItem) -> Result<Entry> {
    match item {
        HierarchyPageItem::Point(entry) => Ok(*entry),
        HierarchyPageItem::Child(child) => Ok(Entry {
            key: child.key,
            offset: child.offset,
            byte_size: i32::try_from(child.byte_size).map_err(|_| {
                Error::InvalidInput("child hierarchy page exceeds COPC i32 byte size".into())
            })?,
            point_count: -1,
        }),
    }
}

fn key_contains(ancestor: VoxelKey, key: VoxelKey) -> bool {
    if key.level < ancestor.level {
        return false;
    }
    let shift = (key.level - ancestor.level) as u32;
    (key.x >> shift) == ancestor.x
        && (key.y >> shift) == ancestor.y
        && (key.z >> shift) == ancestor.z
}

struct LodIndex {
    nodes: Vec<LodNodeRange>,
    order_path: TempPath,
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct LodNodeRange {
    key: VoxelKey,
    start: u64,
    count: usize,
}

struct IndexRun {
    path: TempPath,
    start: u64,
    count: usize,
}

fn build_lod_index<S: CopcPointSource>(
    source: &S,
    center: (f64, f64, f64),
    halfsize: f64,
    params: &CopcWriterParams,
    cancel: &dyn CancelCheck,
) -> Result<LodIndex> {
    cancel.check()?;
    let total_points = u32::try_from(source.len()).map_err(|_| {
        Error::InvalidInput("COPC writer supports at most u32::MAX points per file".into())
    })?;
    let max_points_per_node = params.max_points_per_node.max(1) as usize;
    let max_depth = params.max_depth.min(30);
    let root_run = write_root_index_run(total_points, cancel)?;
    let mut order_file = new_index_tempfile("order")?;
    let mut order_offset = 0;
    let mut nodes = Vec::new();
    {
        let mut order_writer = BufWriter::new(order_file.as_file_mut());
        let mut builder = LodIndexBuilder {
            source,
            max_points_per_node,
            max_depth,
            cancel,
            order_writer: &mut order_writer,
            order_offset: &mut order_offset,
            nodes: &mut nodes,
        };
        builder.assign(VoxelKey::root(), root_run, Bounds::cube(center, halfsize))?;
        order_writer
            .flush()
            .map_err(|e| Error::io("flush LOD index order", e))?;
    }
    nodes.sort_by_key(|node| node.key);
    Ok(LodIndex {
        nodes,
        order_path: order_file.into_temp_path(),
    })
}

struct LodIndexBuilder<'a, S: CopcPointSource, W: Write> {
    source: &'a S,
    max_points_per_node: usize,
    max_depth: u32,
    cancel: &'a dyn CancelCheck,
    order_writer: &'a mut W,
    order_offset: &'a mut u64,
    nodes: &'a mut Vec<LodNodeRange>,
}

impl<S: CopcPointSource, W: Write> LodIndexBuilder<'_, S, W> {
    fn assign(&mut self, key: VoxelKey, run: IndexRun, bounds: Bounds) -> Result<()> {
        self.cancel.check()?;
        if run.count == 0 {
            return Ok(());
        }
        if run.count <= self.max_points_per_node || key.level as u32 >= self.max_depth {
            let start = *self.order_offset;
            append_index_run_to_order(&run, self.order_writer, self.order_offset, self.cancel)?;
            self.nodes.push(LodNodeRange {
                key,
                start,
                count: run.count,
            });
            return Ok(());
        }

        let mut children = partition_index_run(self.source, &run, bounds, self.cancel)?;
        let start = *self.order_offset;
        let selected_counts = append_lod_selection_to_order(
            &children,
            self.max_points_per_node,
            self.order_writer,
            self.order_offset,
            self.cancel,
        )?;
        let selected_total = selected_counts.iter().sum();
        self.nodes.push(LodNodeRange {
            key,
            start,
            count: selected_total,
        });

        for (octant, child) in children.iter_mut().enumerate() {
            let Some(mut child_run) = child.take() else {
                continue;
            };
            let selected = selected_counts[octant];
            if selected >= child_run.count {
                continue;
            }
            child_run.start += selected as u64 * INDEX_RECORD_BYTES;
            child_run.count -= selected;
            self.assign(
                key.child(octant as u8),
                child_run,
                bounds.octant(octant as u8),
            )?;
        }
        Ok(())
    }
}

fn write_root_index_run(total_points: u32, cancel: &dyn CancelCheck) -> Result<IndexRun> {
    let mut writer = BufWriter::new(new_index_tempfile("root")?);
    for index in 0..total_points {
        if index as usize % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        writer
            .write_u32::<LittleEndian>(index)
            .map_err(|e| Error::io("write root LOD index", e))?;
    }
    let file = writer
        .into_inner()
        .map_err(|e| Error::io("flush root LOD index", e.into_error()))?;
    Ok(IndexRun {
        path: file.into_temp_path(),
        start: 0,
        count: total_points as usize,
    })
}

fn partition_index_run<S: CopcPointSource>(
    source: &S,
    run: &IndexRun,
    bounds: Bounds,
    cancel: &dyn CancelCheck,
) -> Result<[Option<IndexRun>; 8]> {
    let mut reader = open_index_run(run)?;
    let mut writers: [Option<BufWriter<NamedTempFile>>; 8] = std::array::from_fn(|_| None);
    let mut counts = [0usize; 8];
    for read_index in 0..run.count {
        if read_index % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        let index = reader
            .read_u32::<LittleEndian>()
            .map_err(|e| Error::io("read LOD partition index", e))?;
        let (x, y, z) = source.xyz(index as usize);
        let octant = child_octant(bounds, x, y, z);
        if writers[octant].is_none() {
            writers[octant] = Some(BufWriter::new(new_index_tempfile("partition")?));
        }
        writers[octant]
            .as_mut()
            .expect("partition writer exists")
            .write_u32::<LittleEndian>(index)
            .map_err(|e| Error::io("write LOD partition index", e))?;
        counts[octant] += 1;
    }

    let mut children: [Option<IndexRun>; 8] = std::array::from_fn(|_| None);
    for octant in 0..8 {
        let Some(writer) = writers[octant].take() else {
            continue;
        };
        let file = writer
            .into_inner()
            .map_err(|e| Error::io("flush LOD partition index", e.into_error()))?;
        children[octant] = Some(IndexRun {
            path: file.into_temp_path(),
            start: 0,
            count: counts[octant],
        });
    }
    Ok(children)
}

fn append_lod_selection_to_order<W: Write>(
    children: &[Option<IndexRun>; 8],
    max_points_per_node: usize,
    order_writer: &mut W,
    order_offset: &mut u64,
    cancel: &dyn CancelCheck,
) -> Result<[usize; 8]> {
    let mut readers: [Option<BufReader<File>>; 8] = std::array::from_fn(|_| None);
    for octant in 0..8 {
        if let Some(child) = &children[octant] {
            readers[octant] = Some(open_index_run(child)?);
        }
    }

    let mut selected_counts = [0usize; 8];
    let mut selected_total = 0usize;
    while selected_total < max_points_per_node {
        cancel.check()?;
        let mut progressed = false;
        for octant in 0..8 {
            let Some(child) = &children[octant] else {
                continue;
            };
            if selected_counts[octant] >= child.count {
                continue;
            }
            let index = readers[octant]
                .as_mut()
                .expect("partition reader exists")
                .read_u32::<LittleEndian>()
                .map_err(|e| Error::io("read selected LOD index", e))?;
            append_index_to_order(order_writer, order_offset, index)?;
            selected_counts[octant] += 1;
            selected_total += 1;
            progressed = true;
            if selected_total == max_points_per_node {
                break;
            }
        }
        if !progressed {
            break;
        }
    }
    Ok(selected_counts)
}

fn append_index_run_to_order<W: Write>(
    run: &IndexRun,
    order_writer: &mut W,
    order_offset: &mut u64,
    cancel: &dyn CancelCheck,
) -> Result<()> {
    let mut reader = open_index_run(run)?;
    for read_index in 0..run.count {
        if read_index % CANCEL_POLL_STRIDE == 0 {
            cancel.check()?;
        }
        let index = reader
            .read_u32::<LittleEndian>()
            .map_err(|e| Error::io("read LOD index", e))?;
        append_index_to_order(order_writer, order_offset, index)?;
    }
    Ok(())
}

fn append_index_to_order<W: Write>(
    order_writer: &mut W,
    order_offset: &mut u64,
    index: u32,
) -> Result<()> {
    order_writer
        .write_u32::<LittleEndian>(index)
        .map_err(|e| Error::io("write LOD index order", e))?;
    *order_offset = order_offset
        .checked_add(INDEX_RECORD_BYTES)
        .ok_or_else(|| Error::InvalidInput("LOD index order exceeds u64 range".into()))?;
    Ok(())
}

fn open_index_run(run: &IndexRun) -> Result<BufReader<File>> {
    let path: &Path = run.path.as_ref();
    let mut file = File::open(path).map_err(|e| Error::io("open LOD index", e))?;
    file.seek(SeekFrom::Start(run.start))
        .map_err(|e| Error::io("seek LOD index", e))?;
    Ok(BufReader::new(file))
}

fn new_index_tempfile(label: &str) -> Result<NamedTempFile> {
    let prefix = format!(".copc-writer-{label}.");
    tempfile::Builder::new()
        .prefix(&prefix)
        .suffix(".idx")
        .tempfile()
        .map_err(|e| Error::io("create LOD index file", e))
}

fn child_octant(bounds: Bounds, x: f64, y: f64, z: f64) -> usize {
    let center = bounds.center();
    usize::from(x >= center.0)
        | (usize::from(y >= center.1) << 1)
        | (usize::from(z >= center.2) << 2)
}

fn cube_from_bounds(bounds: &Bounds) -> ((f64, f64, f64), f64) {
    let dx = bounds.max.0 - bounds.min.0;
    let dy = bounds.max.1 - bounds.min.1;
    let dz = bounds.max.2 - bounds.min.2;
    let center = (
        bounds.min.0 + dx * 0.5,
        bounds.min.1 + dy * 0.5,
        bounds.min.2 + dz * 0.5,
    );
    let halfsize = (dx.max(dy).max(dz) * 0.5).max(1e-6);
    (center, halfsize)
}

struct LasHeader {
    point_data_format: u8,
    point_record_length: u16,
    offset_to_point_data: u32,
    number_of_vlrs: u32,
    file_source_id: u16,
    global_encoding: u16,
    guid: [u8; 16],
    system_identifier: String,
    generating_software: String,
    creation_day_of_year: u16,
    creation_year: u16,
    scale: (f64, f64, f64),
    offset: (f64, f64, f64),
    bounds: Bounds,
    legacy_point_count: u32,
    total_point_count: u64,
    offset_to_first_evlr: u64,
    number_of_evlrs: u32,
    extended_return_counts: [u64; 15],
}

impl LasHeader {
    fn write<W: Write>(&self, writer: &mut W) -> Result<()> {
        writer
            .write_all(b"LASF")
            .map_err(|e| Error::io("write LAS signature", e))?;
        writer
            .write_u16::<LittleEndian>(self.file_source_id)
            .map_err(|e| Error::io("write file source id", e))?;
        writer
            .write_u16::<LittleEndian>(self.global_encoding)
            .map_err(|e| Error::io("write global encoding", e))?;
        writer
            .write_all(&self.guid)
            .map_err(|e| Error::io("write GUID", e))?;
        writer
            .write_u8(1)
            .map_err(|e| Error::io("write version major", e))?;
        writer
            .write_u8(4)
            .map_err(|e| Error::io("write version minor", e))?;
        writer
            .write_all(&pad(self.system_identifier.as_bytes(), 32))
            .map_err(|e| Error::io("write system id", e))?;
        writer
            .write_all(&pad(self.generating_software.as_bytes(), 32))
            .map_err(|e| Error::io("write generating software", e))?;
        writer
            .write_u16::<LittleEndian>(self.creation_day_of_year)
            .map_err(|e| Error::io("write creation day", e))?;
        writer
            .write_u16::<LittleEndian>(self.creation_year)
            .map_err(|e| Error::io("write creation year", e))?;
        writer
            .write_u16::<LittleEndian>(375)
            .map_err(|e| Error::io("write header size", e))?;
        writer
            .write_u32::<LittleEndian>(self.offset_to_point_data)
            .map_err(|e| Error::io("write point data offset", e))?;
        writer
            .write_u32::<LittleEndian>(self.number_of_vlrs)
            .map_err(|e| Error::io("write VLR count", e))?;
        writer
            .write_u8(self.point_data_format)
            .map_err(|e| Error::io("write point format", e))?;
        writer
            .write_u16::<LittleEndian>(self.point_record_length)
            .map_err(|e| Error::io("write point record length", e))?;
        writer
            .write_u32::<LittleEndian>(self.legacy_point_count)
            .map_err(|e| Error::io("write legacy point count", e))?;
        for _ in 0..5 {
            writer
                .write_u32::<LittleEndian>(0)
                .map_err(|e| Error::io("write legacy returns", e))?;
        }
        writer
            .write_f64::<LittleEndian>(self.scale.0)
            .map_err(|e| Error::io("write x scale", e))?;
        writer
            .write_f64::<LittleEndian>(self.scale.1)
            .map_err(|e| Error::io("write y scale", e))?;
        writer
            .write_f64::<LittleEndian>(self.scale.2)
            .map_err(|e| Error::io("write z scale", e))?;
        writer
            .write_f64::<LittleEndian>(self.offset.0)
            .map_err(|e| Error::io("write x offset", e))?;
        writer
            .write_f64::<LittleEndian>(self.offset.1)
            .map_err(|e| Error::io("write y offset", e))?;
        writer
            .write_f64::<LittleEndian>(self.offset.2)
            .map_err(|e| Error::io("write z offset", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.max.0)
            .map_err(|e| Error::io("write max x", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.min.0)
            .map_err(|e| Error::io("write min x", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.max.1)
            .map_err(|e| Error::io("write max y", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.min.1)
            .map_err(|e| Error::io("write min y", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.max.2)
            .map_err(|e| Error::io("write max z", e))?;
        writer
            .write_f64::<LittleEndian>(self.bounds.min.2)
            .map_err(|e| Error::io("write min z", e))?;
        writer
            .write_u64::<LittleEndian>(0)
            .map_err(|e| Error::io("write waveform packet start", e))?;
        writer
            .write_u64::<LittleEndian>(self.offset_to_first_evlr)
            .map_err(|e| Error::io("write first EVLR offset", e))?;
        writer
            .write_u32::<LittleEndian>(self.number_of_evlrs)
            .map_err(|e| Error::io("write EVLR count", e))?;
        writer
            .write_u64::<LittleEndian>(self.total_point_count)
            .map_err(|e| Error::io("write total point count", e))?;
        for count in self.extended_return_counts {
            writer
                .write_u64::<LittleEndian>(count)
                .map_err(|e| Error::io("write extended returns", e))?;
        }
        Ok(())
    }
}

fn pad(value: &[u8], len: usize) -> Vec<u8> {
    let mut out = Vec::with_capacity(len);
    let take = value.len().min(len);
    out.extend_from_slice(&value[..take]);
    out.resize(len, 0);
    out
}

fn write_vlr_header<W: Write>(
    writer: &mut W,
    user_id: &str,
    record_id: u16,
    body_size: u16,
    description: &str,
) -> Result<()> {
    writer
        .write_u16::<LittleEndian>(0)
        .map_err(|e| Error::io("write VLR reserved", e))?;
    writer
        .write_all(&pad(user_id.as_bytes(), 16))
        .map_err(|e| Error::io("write VLR user id", e))?;
    writer
        .write_u16::<LittleEndian>(record_id)
        .map_err(|e| Error::io("write VLR record id", e))?;
    writer
        .write_u16::<LittleEndian>(body_size)
        .map_err(|e| Error::io("write VLR body size", e))?;
    writer
        .write_all(&pad(description.as_bytes(), 32))
        .map_err(|e| Error::io("write VLR description", e))?;
    Ok(())
}

fn write_evlr_header<W: Write>(
    writer: &mut W,
    user_id: &str,
    record_id: u16,
    body_size: u64,
    description: &str,
) -> Result<()> {
    writer
        .write_u16::<LittleEndian>(0)
        .map_err(|e| Error::io("write EVLR reserved", e))?;
    writer
        .write_all(&pad(user_id.as_bytes(), 16))
        .map_err(|e| Error::io("write EVLR user id", e))?;
    writer
        .write_u16::<LittleEndian>(record_id)
        .map_err(|e| Error::io("write EVLR record id", e))?;
    writer
        .write_u64::<LittleEndian>(body_size)
        .map_err(|e| Error::io("write EVLR body size", e))?;
    writer
        .write_all(&pad(description.as_bytes(), 32))
        .map_err(|e| Error::io("write EVLR description", e))?;
    Ok(())
}

fn encode_point_record(
    buf: &mut [u8],
    fields: &CopcPointFields,
    scale: (f64, f64, f64),
    offset: (f64, f64, f64),
    point_index: usize,
    format: &LasFormat,
) -> Result<()> {
    let mut cursor = Cursor::new(buf);
    let (ix, iy, iz) = quantize_xyz(point_index, fields.x, fields.y, fields.z, scale, offset)?;
    let flags =
        fields.synthetic | (fields.key_point << 1) | (fields.withheld << 2) | (fields.overlap << 3);
    let point = raw::Point {
        x: ix,
        y: iy,
        z: iz,
        intensity: fields.intensity,
        flags: raw::point::Flags::ThreeByte(
            fields.return_number | (fields.number_of_returns << 4),
            flags
                | (fields.scan_channel << 4)
                | (fields.scan_direction_flag << 6)
                | (fields.edge_of_flight_line << 7),
            fields.classification,
        ),
        scan_angle: raw::point::ScanAngle::from(fields.scan_angle),
        user_data: fields.user_data,
        point_source_id: fields.point_source_id,
        gps_time: Some(fields.gps_time),
        color: format
            .has_color
            .then_some(Color::new(fields.red, fields.green, fields.blue)),
        waveform: None,
        nir: None,
        extra_bytes: Vec::new(),
    };
    point
        .write_to(&mut cursor, format)
        .map_err(|e| Error::Las(format!("write point record: {e}")))?;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    struct VecSource {
        points: Vec<CopcPointFields>,
    }

    impl CopcPointSource for VecSource {
        fn len(&self) -> usize {
            self.points.len()
        }

        fn xyz(&self, index: usize) -> (f64, f64, f64) {
            let point = self.points[index];
            (point.x, point.y, point.z)
        }

        fn fields(&self, index: usize) -> Result<CopcPointFields> {
            Ok(self.points[index])
        }
    }

    #[test]
    fn spooled_lod_index_covers_each_point_once() {
        let points = (0..257)
            .map(|i| CopcPointFields {
                x: f64::from((i * 37) % 101),
                y: f64::from((i * 53) % 103),
                z: f64::from((i * 71) % 107),
                intensity: 0,
                return_number: 1,
                number_of_returns: 1,
                synthetic: 0,
                key_point: 0,
                withheld: 0,
                overlap: 0,
                scan_channel: 0,
                scan_direction_flag: 0,
                edge_of_flight_line: 0,
                classification: 0,
                user_data: 0,
                scan_angle: 0.0,
                point_source_id: 0,
                gps_time: f64::from(i),
                red: 0,
                green: 0,
                blue: 0,
            })
            .collect();
        let source = VecSource { points };
        let bounds = source_bounds(&source);
        let (center, halfsize) = cube_from_bounds(&bounds);
        let params = CopcWriterParams {
            max_points_per_node: 7,
            max_depth: 5,
        };

        let spooled = build_lod_index(&source, center, halfsize, &params, &NeverCancel).unwrap();
        let ranges = read_lod_index(&spooled).unwrap();

        let mut seen = vec![false; source.len()];
        let mut total = 0usize;
        for (key, indices) in ranges {
            if key.level < params.max_depth as i32 {
                assert!(indices.len() <= params.max_points_per_node as usize);
            }
            for index in indices {
                let seen = &mut seen[index as usize];
                assert!(!*seen, "point index {index} was assigned more than once");
                *seen = true;
                total += 1;
            }
        }
        assert_eq!(source.len(), total);
        assert!(seen.into_iter().all(|value| value));
    }

    #[test]
    fn hierarchy_plan_splits_large_root_page() {
        let mut entries = vec![Entry {
            key: VoxelKey::root(),
            offset: 1,
            byte_size: 1,
            point_count: 1,
        }];
        let mut offset = 2;
        for z in 0..16 {
            for y in 0..16 {
                for x in 0..16 {
                    entries.push(Entry {
                        key: VoxelKey { level: 4, x, y, z },
                        offset,
                        byte_size: 1,
                        point_count: 1,
                    });
                    offset += 1;
                }
            }
        }
        entries.sort_by_key(|entry| entry.key);

        let mut plan = plan_hierarchy_pages(&entries, VoxelKey::root()).unwrap();
        let start = 1024;
        let end = assign_hierarchy_page_offsets(&mut plan, start).unwrap();

        assert!(plan.byte_size < hierarchy_page_byte_size(entries.len()).unwrap());
        assert!(plan
            .items
            .iter()
            .any(|item| matches!(item, HierarchyPageItem::Child(_))));

        let mut out = Cursor::new(vec![0; start as usize]);
        out.seek(SeekFrom::Start(start)).unwrap();
        write_hierarchy_page_tree(&mut out, &plan).unwrap();
        assert_eq!(end, out.get_ref().len() as u64);
    }

    fn source_bounds(source: &VecSource) -> Bounds {
        source.points.iter().fold(
            Bounds::point(source.points[0].x, source.points[0].y, source.points[0].z),
            |mut bounds, point| {
                bounds.extend(point.x, point.y, point.z);
                bounds
            },
        )
    }

    fn read_lod_index(index: &LodIndex) -> Result<Vec<(VoxelKey, Vec<u32>)>> {
        let path: &Path = index.order_path.as_ref();
        let mut reader =
            BufReader::new(File::open(path).map_err(|e| Error::io("open LOD order", e))?);
        let mut out = Vec::new();
        for node in &index.nodes {
            reader
                .seek(SeekFrom::Start(node.start))
                .map_err(|e| Error::io("seek LOD order", e))?;
            let mut indices = Vec::with_capacity(node.count);
            for _ in 0..node.count {
                indices.push(
                    reader
                        .read_u32::<LittleEndian>()
                        .map_err(|e| Error::io("read LOD order", e))?,
                );
            }
            out.push((node.key, indices));
        }
        Ok(out)
    }
}