From d56458cf0b958af2d86c391d1eea140e43d4d5b3 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 10 Jan 2018 17:32:36 +0100 Subject: [PATCH 1/3] [WIP] First DataStreams version without LasPoint initialization. --- REQUIRE | 2 + src/LasIO.jl | 4 ++ src/datastreams.jl | 119 +++++++++++++++++++++++++++++++++++++++++++++ src/fileio.jl | 17 ++++++- src/header.jl | 68 +++++++++++++------------- src/point.jl | 76 +++++++++++++++-------------- 6 files changed, 214 insertions(+), 72 deletions(-) create mode 100644 src/datastreams.jl diff --git a/REQUIRE b/REQUIRE index 1dd68ec..8b550ca 100644 --- a/REQUIRE +++ b/REQUIRE @@ -3,3 +3,5 @@ ColorTypes 0.3.3 FixedPointNumbers 0.3.2 FileIO 0.3.0 GeometryTypes 0.4.0 +DataStreams +Parameters \ No newline at end of file diff --git a/src/LasIO.jl b/src/LasIO.jl index 03ff034..5b84e82 100644 --- a/src/LasIO.jl +++ b/src/LasIO.jl @@ -2,10 +2,12 @@ __precompile__() module LasIO +using DataStreams using FileIO using FixedPointNumbers using ColorTypes using GeometryTypes # for conversion +using Parameters export # Types @@ -19,6 +21,7 @@ export # Functions on LasHeader update!, + pointtype, # Functions on LasPoint return_number, @@ -52,6 +55,7 @@ include("point.jl") include("util.jl") include("fileio.jl") include("srs.jl") +include("datastreams.jl") function __init__() # these should eventually go in diff --git a/src/datastreams.jl b/src/datastreams.jl new file mode 100644 index 0000000..33df8da --- /dev/null +++ b/src/datastreams.jl @@ -0,0 +1,119 @@ +######## +# SOURCE +######## + + +mutable struct Source <: Data.Source + schema::Data.Schema + header::LasHeader + io::IO + fullpath::String + datapos::Int +end + +dict_of_struct(T) = Dict((String(fieldname(typeof(T), i)), getfield(T, i)) for i = 1:nfields(T)) + +function Source(f::AbstractString) + # s = is_windows() ? open(f) : IOStream(Mmap.mmap(f)) + io = open(f) + + skiplasf(io) + header = read(io, LasHeader) + + n = header.records_count + pointtype = pointformat(header) + + # Schema + columns = Vector{String}(fieldnames(pointtype)) + types = Vector{Type}([fieldtype(pointtype, i) for i in 1:nfields(pointtype)]) + sch = Data.Schema(types, columns, n, dict_of_struct(header)) + + return Source(sch, header, io, f, position(io)) +end + +Data.reset!(s::LasIO.Source) = (seek(s.io, s.datapos); return nothing) +Data.schema(source::LasIO.Source) = source.schema +Data.accesspattern(::Type{LasIO.Source}) = Data.Sequential +Data.isdone(source::LasIO.Source, row::Int, col::Int) = eof(source.io) || (row, col) > Data.size(Data.schema(source)) +Data.isdone(source::LasIO.Source, row, col, rows, cols) = eof(source.io) || row > rows || col > cols +Data.streamtype(::Type{<:LasIO.Source}, ::Type{Data.Field}) = true + +# Data.streamfrom(source::LasIO.Source, st::Type{Data.Row}, t::Type{T}, row::Int) where {T} = read(source.io, pointformat(source.header)) +function Data.streamfrom(source::LasIO.Source, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T} + # p = source.pointdata[row] + # v = getfield(p, fieldname(typeof(p), col)) + read(source.io, T) +end + +###### +# SINK +###### + +mutable struct Sink{T} <: Data.Sink where T <: LasPoint + stream::IO + header::LasHeader + # pointdata::Vector{T} + pointformat::Type{T} +end + +# setup header and empty pointvector +function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::AbstractString) + s = open(fullpath, "w") + + # determine point version and derivatives + pointtype = gettype(Data.header(sch)) + data_format_id = pointformat(pointtype) + data_record_length = packed_size(pointtype) + n = Data.size(sch, 1) # rows + + # create header + header = LasHeader(data_format_id=data_format_id, data_record_length=data_record_length, records_count=n) + + # write header + write(s, magic(format"LAS")) + write(s, header) + println("Stream now at $(position(s))") + + # create empty pointdata + # pointdata = Vector{pointtype}(n) + + # return stream, position is correct for writing points + return Sink(s, header, pointtype) +end + +# Update existing Sink +# function Sink(sink::LasIO.Sink, sch::Data.Schema, S::Type{StreamType}, append::Bool; reference::Vector{UInt8}=UInt8[]) + # return Sink +# end + +"""Determine LAS versions based on specific columns.""" +function gettype(columns::Vector{String}) + # LAS versions only differ in gps_time and rgb color information + has_gps, has_color = false, false + ("gps_time" in columns) && (has_gps = true) + ("red" in columns || "green" in columns || "blue" in columns) && (has_color = true) + + has_gps && has_color && return LasPoint3 + has_color && return LasPoint2 + has_gps && return LasPoint1 + return LasPoint0 +end + +Data.streamtypes(::Type{LasIO.Sink}) = [Data.Field, Data.Row] +Data.cleanup!(sink::LasIO.Sink) = nothing +Data.weakrefstrings(::Type{LasIO.Sink}) = false + +# actually write points to our pointvector +function Data.streamto!(sink::LasIO.Sink, S::Type{Data.Field}, val, row, col) + # TODO(evetion) check if stream position is at row*col + # write points + write(sink.stream, val) +end + +# save file +function Data.close!(sink::LasIO.Sink) + close(sink.stream) + # TODO(evetion) check number of points, extents etc + # change header if possible + return sink +end diff --git a/src/fileio.jl b/src/fileio.jl index a540ce2..a30d4e5 100644 --- a/src/fileio.jl +++ b/src/fileio.jl @@ -14,6 +14,20 @@ function pointformat(header::LasHeader) end end +function pointformat(t::Type{T}) where T <: LasPoint + if t == LasPoint0 + return 0x00 + elseif t == LasPoint1 + return 0x01 + elseif t == LasPoint2 + return 0x02 + elseif t == LasPoint3 + return 0x03 + else + error("unsupported point format $t") + end +end + # skip the LAS file's magic four bytes, "LASF" skiplasf(s::Union{Stream{format"LAS"}, Stream{format"LAZ"}, IO}) = read(s, UInt32) @@ -45,8 +59,7 @@ end function read_header(f::AbstractString) open(f) do s - skiplasf(s) - read(s, LasHeader) + read_header(s::IO) end end diff --git a/src/header.jl b/src/header.jl index c5774f9..a210ca3 100644 --- a/src/header.jl +++ b/src/header.jl @@ -14,40 +14,40 @@ Backward compatibility with LAS 1.1 – LAS 1.3 when payloads consist of only le content =# -mutable struct LasHeader - file_source_id::UInt16 - global_encoding::UInt16 - guid_1::UInt32 - guid_2::UInt16 - guid_3::UInt16 - guid_4::AbstractString - version_major::UInt8 - version_minor::UInt8 - system_id::AbstractString - software_id::AbstractString - creation_doy::UInt16 - creation_year::UInt16 - header_size::UInt16 - data_offset::UInt32 - n_vlr::UInt32 - data_format_id::UInt8 - data_record_length::UInt16 - records_count::UInt32 - point_return_count::Vector{UInt32} - x_scale::Float64 - y_scale::Float64 - z_scale::Float64 - x_offset::Float64 - y_offset::Float64 - z_offset::Float64 - x_max::Float64 - x_min::Float64 - y_max::Float64 - y_min::Float64 - z_max::Float64 - z_min::Float64 - variable_length_records::Vector{LasVariableLengthRecord} - user_defined_bytes::Vector{UInt8} +@with_kw mutable struct LasHeader + file_source_id::UInt16 = UInt16(0) # no known flightline + global_encoding::UInt16 = UInt16(0) + guid_1::UInt32 = UInt32(0) # project IDs + guid_2::UInt16 = UInt16(0) + guid_3::UInt16 = UInt16(0) + guid_4::AbstractString = rpad("", 8) + version_major::UInt8 = UInt8(1) + version_minor::UInt8 = UInt8(0) + system_id::AbstractString = rpad("LasIO.jl datastream", 32) # small char array would be nicer + software_id::AbstractString = rpad("LasIO.jl", 32) + creation_doy::UInt16 = UInt16(Dates.dayofyear(now())) + creation_year::UInt16 = UInt16(Dates.year(now())) + header_size::UInt16 = UInt16(227) + data_offset::UInt32 = UInt32(227) + n_vlr::UInt32 = UInt32(0) + data_format_id::UInt8 = UInt8(pointformat(LasPoint0)) + data_record_length::UInt16 = UInt16(packed_size(LasPoint0)) + records_count::UInt32 = UInt32(0) + point_return_count::Vector{UInt32} = Vector{UInt32}([0,0,0,0,0]) + x_scale::Float64 = 0.01 + y_scale::Float64 = 0.01 + z_scale::Float64 = 0.01 + x_offset::Float64 = 0.0 + y_offset::Float64 = 0.0 + z_offset::Float64 = 0.0 + x_max::Float64 = 0.0 + x_min::Float64 = 0.0 + y_max::Float64 = 0.0 + y_min::Float64 = 0.0 + z_max::Float64 = 0.0 + z_min::Float64 = 0.0 + variable_length_records::Vector{LasVariableLengthRecord} = Vector{LasVariableLengthRecord}() + user_defined_bytes::Vector{UInt8} = Vector{UInt8}() end function Base.show(io::IO, header::LasHeader) diff --git a/src/point.jl b/src/point.jl index 2db6cac..7606fdb 100644 --- a/src/point.jl +++ b/src/point.jl @@ -8,64 +8,68 @@ function Base.show(io::IO, pointdata::Vector{T}) where T <: LasPoint end "ASPRS LAS point data record format 0" -struct LasPoint0 <: LasPoint +@with_kw struct LasPoint0 <: LasPoint x::Int32 y::Int32 z::Int32 - intensity::UInt16 - flag_byte::UInt8 - raw_classification::UInt8 - scan_angle::Int8 - user_data::UInt8 - pt_src_id::UInt16 + intensity::UInt16 = UInt16(0) + flag_byte::UInt8 = 0x0 + raw_classification::UInt8 = 0x0 + scan_angle::Int8 = Int8(0) + user_data::UInt8 = 0x0 + pt_src_id::UInt16 = UInt16(0) end +packed_size(::Type{LasPoint0}) = 20 "ASPRS LAS point data record format 1" -struct LasPoint1 <: LasPoint +@with_kw struct LasPoint1 <: LasPoint x::Int32 y::Int32 z::Int32 - intensity::UInt16 - flag_byte::UInt8 - raw_classification::UInt8 - scan_angle::Int8 - user_data::UInt8 - pt_src_id::UInt16 - gps_time::Float64 + intensity::UInt16 = UInt16(0) + flag_byte::UInt8 = 0x0 + raw_classification::UInt8 = 0x0 + scan_angle::Int8 = Int8(0) + user_data::UInt8 = 0x0 + pt_src_id::UInt16 = UInt16(0) + gps_time::Float64 = gps_time(now()) end +packed_size(::Type{LasPoint1}) = 28 "ASPRS LAS point data record format 2" -struct LasPoint2 <: LasPoint +@with_kw struct LasPoint2 <: LasPoint x::Int32 y::Int32 z::Int32 - intensity::UInt16 - flag_byte::UInt8 - raw_classification::UInt8 - scan_angle::Int8 - user_data::UInt8 - pt_src_id::UInt16 - red::N0f16 - green::N0f16 - blue::N0f16 + intensity::UInt16 = UInt16(0) + flag_byte::UInt8 = 0x0 + raw_classification::UInt8 = 0x0 + scan_angle::Int8 = Int8(0) + user_data::UInt8 = 0x0 + pt_src_id::UInt16 = UInt16(0) + red::N0f16 = N0f16(0) + green::N0f16 = N0f16(0) + blue::N0f16 = N0f16(0) end +packed_size(::Type{LasPoint2}) = 26 "ASPRS LAS point data record format 3" -struct LasPoint3 <: LasPoint +@with_kw struct LasPoint3 <: LasPoint x::Int32 y::Int32 z::Int32 - intensity::UInt16 - flag_byte::UInt8 - raw_classification::UInt8 - scan_angle::Int8 - user_data::UInt8 - pt_src_id::UInt16 - gps_time::Float64 - red::N0f16 - green::N0f16 - blue::N0f16 + intensity::UInt16 = UInt16(0) + flag_byte::UInt8 = 0x0 + raw_classification::UInt8 = 0x0 + scan_angle::Int8 = Int8(0) + user_data::UInt8 = 0x0 + pt_src_id::UInt16 = UInt16(0) + gps_time::Float64 = gps_time(now()) + red::N0f16 = N0f16(0) + green::N0f16 = N0f16(0) + blue::N0f16 = N0f16(0) end +packed_size(::Type{LasPoint3}) = 34 # for convenience in function signatures const LasPointColor = Union{LasPoint2,LasPoint3} From 66f82861fa163765694657b87d99c7bc637bfaf2 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 10 Jan 2018 20:11:33 +0100 Subject: [PATCH 2/3] Header is now updated by streaming points to Sink. --- src/datastreams.jl | 54 ++++++++++++++++++++++++++++++++++++++++------ src/fileio.jl | 2 +- 2 files changed, 48 insertions(+), 8 deletions(-) diff --git a/src/datastreams.jl b/src/datastreams.jl index 33df8da..9648d35 100644 --- a/src/datastreams.jl +++ b/src/datastreams.jl @@ -14,8 +14,7 @@ end dict_of_struct(T) = Dict((String(fieldname(typeof(T), i)), getfield(T, i)) for i = 1:nfields(T)) function Source(f::AbstractString) - # s = is_windows() ? open(f) : IOStream(Mmap.mmap(f)) - io = open(f) + io = is_windows() ? open(f) : IOBuffer(Mmap.mmap(f)) skiplasf(io) header = read(io, LasHeader) @@ -52,8 +51,9 @@ end mutable struct Sink{T} <: Data.Sink where T <: LasPoint stream::IO header::LasHeader - # pointdata::Vector{T} pointformat::Type{T} + bbox::Vector{Float64} + returncount::Vector{UInt32} end # setup header and empty pointvector @@ -77,8 +77,14 @@ function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::Abs # create empty pointdata # pointdata = Vector{pointtype}(n) + # empty bbox (x, x, y, y, z, z) + bbox = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + + # empty return count + return_count = Array{UInt32}([0, 0, 0, 0, 0]) + # return stream, position is correct for writing points - return Sink(s, header, pointtype) + return Sink(s, header, pointtype, bbox, return_count) end # Update existing Sink @@ -103,17 +109,51 @@ Data.streamtypes(::Type{LasIO.Sink}) = [Data.Field, Data.Row] Data.cleanup!(sink::LasIO.Sink) = nothing Data.weakrefstrings(::Type{LasIO.Sink}) = false +function update_sink(sink::LasIO.Sink, col::Integer, val::T) where {T} + if col == 1 + (val > sink.bbox[1] || sink.bbox[1] == 0.0) && (setindex!(sink.bbox, val, 1)) # xmax + (val < sink.bbox[2] || sink.bbox[2] == 0.0) && (setindex!(sink.bbox, val, 2)) # xmin + end + if col == 2 + (val > sink.bbox[3] || sink.bbox[3] == 0.0) && (setindex!(sink.bbox, val, 3)) # ymax + (val < sink.bbox[4] || sink.bbox[4] == 0.0) && (setindex!(sink.bbox, val, 4)) # ymin + end + if col == 3 + (val > sink.bbox[5] || sink.bbox[5] == 0.0) && (setindex!(sink.bbox, val, 5)) # zmax + (val < sink.bbox[6] || sink.bbox[6] == 0.0) && (setindex!(sink.bbox, val, 6)) # zmin + end + if col == 5 + return_number = val & 0b00000111 + return_number <= 5 && (sink.returncount[return_number+1] += 1) + end +end + # actually write points to our pointvector function Data.streamto!(sink::LasIO.Sink, S::Type{Data.Field}, val, row, col) # TODO(evetion) check if stream position is at row*col - # write points + update_sink(sink, col, val) write(sink.stream, val) end # save file function Data.close!(sink::LasIO.Sink) + + # update header + # only works when Int32 are passed, already scaled and offset + header = sink.header + header.x_max = muladd(sink.bbox[1], header.x_scale, header.x_offset) + header.x_min = muladd(sink.bbox[2], header.x_scale, header.x_offset) + header.y_max = muladd(sink.bbox[3], header.y_scale, header.y_offset) + header.y_min = muladd(sink.bbox[4], header.y_scale, header.y_offset) + header.z_max = muladd(sink.bbox[5], header.z_scale, header.z_offset) + header.z_min = muladd(sink.bbox[6], header.z_scale, header.z_offset) + header.point_return_count = sink.returncount + + # seek back to beginning and write header + seekstart(sink.stream) + skiplasf(sink.stream) + write(sink.stream, header) + close(sink.stream) - # TODO(evetion) check number of points, extents etc - # change header if possible return sink end diff --git a/src/fileio.jl b/src/fileio.jl index a30d4e5..a34ed5b 100644 --- a/src/fileio.jl +++ b/src/fileio.jl @@ -29,7 +29,7 @@ function pointformat(t::Type{T}) where T <: LasPoint end # skip the LAS file's magic four bytes, "LASF" -skiplasf(s::Union{Stream{format"LAS"}, Stream{format"LAZ"}, IO}) = read(s, UInt32) +skiplasf(s::Union{Stream{format"LAS"}, Stream{format"LAZ"}, IO}) = skip(s, sizeof(UInt32)) function load(f::File{format"LAS"}) open(f) do s From b9ab16f884a78c83711811da0671b5290af89e0c Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Mon, 15 Jan 2018 18:05:45 +0100 Subject: [PATCH 3/3] Coordinate conversion between Float and Int implemented. --- README.md | 48 +++++++++++++++++++++++++++++- src/datastreams.jl | 74 ++++++++++++++++++++++++++++++++-------------- src/util.jl | 6 ++++ 3 files changed, 105 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index d65b831..1f8c56e 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,53 @@ Julia package for reading and writing the LAS lidar format. This is a pure Julia alternative to [LibLAS.jl](https://github.com/visr/LibLAS.jl) or [Laszip.jl](https://github.com/joa-quim/Laszip.jl). Currently only LAS versions 1.1 - 1.3 and point formats 0 - 3 are supported. For LAZ support see below. -If the file fits into memory, it can be loaded using +## Usage + +```julia +julia> using LasIO +julia> s = LasIO.Source("test/srs.las") +LasIO.Source(Data.Schema: +rows: 10 cols: 10 +Columns: + "x" Int32 + "y" Int32 + "z" Int32 + "intensity" UInt16 + "flag_byte" UInt8 + "raw_classification" UInt8 + "scan_angle" Int8 + "user_data" UInt8 + "pt_src_id" UInt16 + "gps_time" Float64, LasHeader with 10 points. +, IOStream(), "test/srs.las", 759) + +julia> using DataStreams, DataFrames +julia> d = Data.stream!(s, DataFrame); +julia> Data.close!(d) +10×10 DataFrames.DataFrame +│ Row │ x │ y │ z │ intensity │ flag_byte │ raw_classification │ scan_angle │ user_data │ pt_src_id │ gps_time │ +├─────┼───────────┼───────────┼────────┼───────────┼───────────┼────────────────────┼────────────┼───────────┼───────────┼───────────┤ +│ 1 │ 2.89814e5 │ 4.32098e6 │ 170.76 │ 0x0104 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 2 │ 2.89815e5 │ 4.32098e6 │ 170.76 │ 0x0118 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 3 │ 2.89815e5 │ 4.32098e6 │ 170.75 │ 0x0118 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 4 │ 2.89816e5 │ 4.32098e6 │ 170.74 │ 0x0118 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 5 │ 2.89816e5 │ 4.32098e6 │ 170.68 │ 0x0104 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 6 │ 2.89817e5 │ 4.32098e6 │ 170.66 │ 0x00f0 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 7 │ 289817.0 │ 4.32098e6 │ 170.63 │ 0x00f0 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 8 │ 2.89818e5 │ 4.32098e6 │ 170.62 │ 0x0118 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 9 │ 289818.0 │ 4.32098e6 │ 170.61 │ 0x0118 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ +│ 10 │ 2.89819e5 │ 4.32098e6 │ 170.58 │ 0x0104 │ 0x30 │ 0x02 │ 0 │ 0x00 │ 0x0000 │ 4.99451e5 │ + +julia> Data.reset!(s) +julia> d = Data.stream!(s, LasIO.Sink, "test_final.las"); +julia> Data.close!(d) +LasIO.Sink{LasIO.LasPoint1}(IOStream(), LasHeader with 10 points. +, LasIO.LasPoint1) +``` + +### Legacy API + +For backwards compatibility, the legacy API is still provided. ```julia using FileIO, LasIO diff --git a/src/datastreams.jl b/src/datastreams.jl index 9648d35..bac7c16 100644 --- a/src/datastreams.jl +++ b/src/datastreams.jl @@ -2,6 +2,8 @@ # SOURCE ######## +const minimum_coordinate = typemin(Int32) +const maximum_coordinate = typemax(Int32) mutable struct Source <: Data.Source schema::Data.Schema @@ -14,6 +16,8 @@ end dict_of_struct(T) = Dict((String(fieldname(typeof(T), i)), getfield(T, i)) for i = 1:nfields(T)) function Source(f::AbstractString) + !isfile(f) && error("Please provide valid path.") + io = is_windows() ? open(f) : IOBuffer(Mmap.mmap(f)) skiplasf(io) @@ -25,6 +29,7 @@ function Source(f::AbstractString) # Schema columns = Vector{String}(fieldnames(pointtype)) types = Vector{Type}([fieldtype(pointtype, i) for i in 1:nfields(pointtype)]) + types[1:3] = [Float64, Float64, Float64] # Convert XYZ coordinates sch = Data.Schema(types, columns, n, dict_of_struct(header)) return Source(sch, header, io, f, position(io)) @@ -39,9 +44,19 @@ Data.streamtype(::Type{<:LasIO.Source}, ::Type{Data.Field}) = true # Data.streamfrom(source::LasIO.Source, st::Type{Data.Row}, t::Type{T}, row::Int) where {T} = read(source.io, pointformat(source.header)) function Data.streamfrom(source::LasIO.Source, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T} - # p = source.pointdata[row] - # v = getfield(p, fieldname(typeof(p), col)) - read(source.io, T) + + # XYZ => scale stored Integers to Floats using header information + if col == 1 + return muladd(read(source.io, Int32), source.header.x_scale, source.header.x_offset) + elseif col == 2 + return muladd(read(source.io, Int32), source.header.y_scale, source.header.y_offset) + elseif col == 3 + return muladd(read(source.io, Int32), source.header.z_scale, source.header.z_offset) + + # Otherwise return read value + else + v = read(source.io, T) + end end ###### @@ -57,7 +72,11 @@ mutable struct Sink{T} <: Data.Sink where T <: LasPoint end # setup header and empty pointvector -function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::AbstractString) +function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::AbstractString, bbox::Vector{Float64}; scale::Real=0.01, epsg=nothing) + + # validate input + length(bbox) != 6 && error("Provide bounding box as (xmin, ymin, zmin, xmax, ymax, zmax)") + s = open(fullpath, "w") # determine point version and derivatives @@ -66,19 +85,21 @@ function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::Abs data_record_length = packed_size(pointtype) n = Data.size(sch, 1) # rows - # create header - header = LasHeader(data_format_id=data_format_id, data_record_length=data_record_length, records_count=n) + # setup and validate scaling based on bbox and scale + x_offset = determine_offset(bbox[1], bbox[4], scale) + y_offset = determine_offset(bbox[2], bbox[5], scale) + z_offset = determine_offset(bbox[3], bbox[6], scale) + + # create header + header = LasHeader(data_format_id=data_format_id, data_record_length=data_record_length, records_count=n, + x_max=bbox[4], x_min=bbox[1], y_max=bbox[5], y_min=bbox[2], z_max=bbox[6], z_min=bbox[3], + x_scale=scale, y_scale=scale, z_scale=scale, + x_offset=x_offset, y_offset=y_offset, z_offset=z_offset) + epsg != nothing && epsg_code!(header, epsg) # write header write(s, magic(format"LAS")) write(s, header) - println("Stream now at $(position(s))") - - # create empty pointdata - # pointdata = Vector{pointtype}(n) - - # empty bbox (x, x, y, y, z, z) - bbox = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # empty return count return_count = Array{UInt32}([0, 0, 0, 0, 0]) @@ -124,7 +145,7 @@ function update_sink(sink::LasIO.Sink, col::Integer, val::T) where {T} end if col == 5 return_number = val & 0b00000111 - return_number <= 5 && (sink.returncount[return_number+1] += 1) + return_number < 5 && (sink.returncount[return_number+1] += 1) end end @@ -132,21 +153,30 @@ end function Data.streamto!(sink::LasIO.Sink, S::Type{Data.Field}, val, row, col) # TODO(evetion) check if stream position is at row*col update_sink(sink, col, val) - write(sink.stream, val) + + # XYZ => scale given Floats to Integers using header information + if col == 1 + write(sink.stream, round(Int32, ((val - sink.header.x_offset) / sink.header.x_scale))) + elseif col == 2 + write(sink.stream, round(Int32, ((val - sink.header.y_offset) / sink.header.y_scale))) + elseif col == 3 + write(sink.stream, round(Int32, ((val - sink.header.z_offset) / sink.header.z_scale))) + else + write(sink.stream, val) + end end # save file function Data.close!(sink::LasIO.Sink) # update header - # only works when Int32 are passed, already scaled and offset header = sink.header - header.x_max = muladd(sink.bbox[1], header.x_scale, header.x_offset) - header.x_min = muladd(sink.bbox[2], header.x_scale, header.x_offset) - header.y_max = muladd(sink.bbox[3], header.y_scale, header.y_offset) - header.y_min = muladd(sink.bbox[4], header.y_scale, header.y_offset) - header.z_max = muladd(sink.bbox[5], header.z_scale, header.z_offset) - header.z_min = muladd(sink.bbox[6], header.z_scale, header.z_offset) + header.x_max = sink.bbox[1] + header.x_min = sink.bbox[2] + header.y_max = sink.bbox[3] + header.y_min = sink.bbox[4] + header.z_max = sink.bbox[5] + header.z_min = sink.bbox[6] header.point_return_count = sink.returncount # seek back to beginning and write header diff --git a/src/util.jl b/src/util.jl index 098dd6d..198c117 100644 --- a/src/util.jl +++ b/src/util.jl @@ -32,3 +32,9 @@ function update!(h::LasHeader, pvec::Vector{T}) where T <: LasPoint h.records_count = length(pvec) nothing end + +"""Determine offset as implemented by LasTools.""" +function determine_offset(min_value::Real, max_value::Real, scale::Real; threshold::Integer=10^7) + s = round(Integer, (min_value + max_value)/scale/2*threshold) # if value rounds to 0, no offset + s * threshold * scale +end