Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add lenunit option to VTK loading #120

Merged
merged 1 commit into from
Sep 24, 2024
Merged

Add lenunit option to VTK loading #120

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/extra/img.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,9 @@ function imgread(fname; lenunit)
values = (; color=vec(data))
# translation followed by rotation is faster
transform = Translate(-dims[1], 0) → Rotate(-π / 2)
# construct the grid
# construct grid
u = lenunit
origin = Point(ntuple(i -> 0.0u, length(dims)))
origin = ntuple(i -> 0.0u, length(dims))
spacing = ntuple(i -> 1.0u, length(dims))
domain = CartesianGrid(dims, origin, spacing) |> transform
georef(values, domain)
Expand Down
20 changes: 12 additions & 8 deletions src/extra/vtkread.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ function vturead(fname; lenunit)
vtk = ReadVTK.VTKFile(fname)

# construct mesh
points = _points(vtk)
points = _points(vtk, lenunit)
connec = _vtuconnec(vtk)
mesh = SimpleMesh(points, connec)

Expand All @@ -59,7 +59,7 @@ function vtpread(fname; lenunit)
vtk = ReadVTK.VTKFile(fname)

# construct mesh
points = _points(vtk)
points = _points(vtk, lenunit)
connec = _vtpconnec(vtk)
mesh = SimpleMesh(points, connec)

Expand All @@ -74,9 +74,11 @@ function vtrread(fname; lenunit)
vtk = ReadVTK.VTKFile(fname)

# construct grid
u = lenunit
coords = ReadVTK.get_coordinates(vtk)
inds = map(!allequal, coords) |> collect
grid = RectilinearGrid(coords[inds]...)
xyz = map(x -> x * u, coords[inds])
grid = RectilinearGrid(xyz...)

# extract data
vtable, etable = _datatables(vtk)
Expand All @@ -89,10 +91,11 @@ function vtsread(fname; lenunit)
vtk = ReadVTK.VTKFile(fname)

# construct grid
u = lenunit
coords = ReadVTK.get_coordinates(vtk)
inds = map(!allequal, coords) |> collect
dims = findall(!, inds) |> Tuple
XYZ = map(A -> dropdims(A; dims), coords[inds])
XYZ = map(A -> dropdims(A; dims) * u, coords[inds])
grid = StructuredGrid(XYZ...)

# extract data
Expand All @@ -106,6 +109,7 @@ function vtiread(fname; lenunit)
vtk = ReadVTK.VTKFile(fname)

# construct grid
u = lenunit
ext = ReadVTK.get_whole_extent(vtk)
# the get_origin and get_spacing functions drop the z dimension if it is empty,
# but the get_whole_extent function does not
Expand All @@ -115,8 +119,8 @@ function vtiread(fname; lenunit)
(ext[2] - ext[1], ext[4] - ext[3], ext[6] - ext[5])
end
inds = findall(!iszero, dims)
origin = ReadVTK.get_origin(vtk) |> Tuple
spacing = ReadVTK.get_spacing(vtk) |> Tuple
origin = Tuple(ReadVTK.get_origin(vtk)) .* u
spacing = Tuple(ReadVTK.get_spacing(vtk)) .* u
grid = CartesianGrid(dims[inds], origin[inds], spacing[inds])

# extract data
Expand All @@ -130,10 +134,10 @@ end
# UTILS
#-------

function _points(vtk)
function _points(vtk, u)
coords = ReadVTK.get_points(vtk)
inds = map(!allequal, eachrow(coords))
[Point(Tuple(c)) for c in eachcol(coords[inds, :])]
[Point(Tuple(c) .* u) for c in eachcol(coords[inds, :])]
end

function _vtuconnec(vtk)
Expand Down
17 changes: 17 additions & 0 deletions test/io/vtk.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,23 @@
vtable = values(gtb, 0)
@test size(eltype(vtable.myVector)) == (2,)
@test eltype(eltype(vtable.myVector)) <: Float32

# custom lenunit
file = ReadVTK.get_example_file("celldata_appended_binary_compressed.vtu", output_directory=savedir)
gtb = GeoIO.load(file, lenunit=cm)
@test unit(Meshes.lentype(crs(gtb.geometry))) == cm
file = joinpath(datadir, "spiral.vtp")
gtb = GeoIO.load(file, lenunit=cm)
@test unit(Meshes.lentype(crs(gtb.geometry))) == cm
file = joinpath(datadir, "rectilinear.vtr")
gtb = GeoIO.load(file, lenunit=cm)
@test unit(Meshes.lentype(crs(gtb.geometry))) == cm
file = joinpath(datadir, "structured.vts")
gtb = GeoIO.load(file, lenunit=cm)
@test unit(Meshes.lentype(crs(gtb.geometry))) == cm
file = joinpath(datadir, "imagedata.vti")
gtb = GeoIO.load(file, lenunit=cm)
@test unit(Meshes.lentype(crs(gtb.geometry))) == cm
end

@testset "save" begin
Expand Down