diff --git a/Project.toml b/Project.toml index 311191b..e5fa1d4 100644 --- a/Project.toml +++ b/Project.toml @@ -32,8 +32,8 @@ PyramidSchemeMakieExt = "Makie" [compat] Aqua = "0.8" ArchGDAL = "0.10.4" -CairoMakie = "0.12,0.13,0.14,0.15" -Colors = "0.12" +CairoMakie = "0.13, 0.14, 0.15" +Colors = "0.12, 0.13" DimensionalData = "0.28, 0.29" DiskArrayEngine = "0.1.2, 0.2" DiskArrayTools = "0.1.10" diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index d0457e8..f709068 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -209,18 +209,31 @@ Fill the pyramids generated from the `data` with the aggregation function `func` `recursive` indicates whether higher tiles are computed from lower tiles or directly from the original data. This is an optimization which for functions like median might lead to misleading results. """ -function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, kwargs...) +function fill_pyramids(data, outputs,func,recursive,spatial_dims;runner=LocalRunner, kwargs...) n_level = length(outputs) pixel_base_size = 2^n_level pyramid_sizes = size.(outputs) tmp_sizes = [ceil(Int,pixel_base_size / 2^i) for i in 1:n_level] - ia = InputArray(data, windows = arraywindows(size(data),pixel_base_size)) + inputwindows = Base.OneTo.(size(data)) + outputwindows = Base.OneTo.(size(data)) + for d in spatial_dims + inputwindows = Base.setindex(inputwindows,RegularWindows(1,size(data,d),window=pixel_base_size),d) + end + + ia = InputArray(data, windows = inputwindows) + - oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = arraywindows(pyramid_sizes[i],tmp_sizes[i])),n_level) + oa = ntuple(n_level) do i + for d in spatial_dims + outputwindows = Base.setindex(outputwindows,RegularWindows(1,pyramid_sizes[i][d],window=tmp_sizes[i]),d) + end + create_outwindows(pyramid_sizes[i],windows = outputwindows) + end func = DiskArrayEngine.create_userfunction(gen_pyr,ntuple(_->eltype(first(outputs)),length(outputs));is_mutating=true,kwargs = (;recursive),args = (func,)) + op = GMDWop((ia,), oa, func) lr = DiskArrayEngine.optimize_loopranges(op,5e8,tol_low=0.2,tol_high=0.05,max_order=2) @@ -317,7 +330,7 @@ output_arrays(pyramid_sizes, T) = [gen_output(T,p) for p in pyramid_sizes] SpatialDim Union of Dimensions which are assumed to be in space and are therefore used in the pyramid building. """ -SpatialDim = Union{DD.Dimensions.XDim, DD.Dimensions.YDim} +SpatialDim = (DD.Dimensions.XDim, DD.Dimensions.YDim) """ buildpyramids(path; resampling_method=mean) @@ -326,7 +339,7 @@ The different scales are written according to the GeoZarr spec and a multiscales The data is aggregated with the specified `resampling_method`. Keyword arguments are forwarded to the `fill_pyramids` function. """ -function buildpyramids(path::AbstractString; resampling_method=mean, recursive=true, runner=LocalRunner, verbose=false) +function buildpyramids(path::AbstractString; resampling_method=mean, recursive=true, runner=LocalRunner, verbose=false, spatial_dims=SpatialDim) if YAB.backendfrompath(path) != YAB.backendlist[:zarr] @show YAB.backendfrompath(path) throw(ArgumentError("$path is not a Zarr dataset therefore we can't build the Pyramids inplace")) @@ -340,17 +353,19 @@ function buildpyramids(path::AbstractString; resampling_method=mean, recursive=t #t = Missing <: eltype(org) ? Union{Missing, tfunc} : tfunc t = Base.infer_return_type(resampling_method, (Matrix{nonmissingtype(eltype(org))},)) - n_level = compute_nlevels(org) - input_axes = filter(x-> x isa SpatialDim, DD.dims(org)) + input_axes = DD.dims(org, spatial_dims) + outarrs = [output_zarr(n, DD.dims(org), t, joinpath(path, string(n)),input_axes) for n in 1:n_level] + if length(input_axes) != 2 throw(ArgumentError("Expected two spatial dimensions got $input_axes")) end verbose && println("Constructing output arrays") - outarrs = [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] verbose && println("Start computation") - fill_pyramids(org, outarrs, resampling_method, recursive;runner) - pyraxs = [agg_axis.(input_axes, 2^n) for n in 1:n_level] + ispatial_dims = DD.dimnum(DD.dims(org),spatial_dims) + fill_pyramids(org, outarrs, resampling_method, recursive, ispatial_dims;runner) + pyraxs_space = [agg_axis.(input_axes, 2^n) for n in 1:n_level] + pyraxs = [DD.setdims(DD.dims(org),p) for p in pyraxs_space] pyrlevels = DD.DimArray.(outarrs, pyraxs) meta = Dict(deepcopy(DD.metadata(org))) push!(meta, "resampling_method" => string(resampling_method)) @@ -370,15 +385,23 @@ end """ output_zarr(n, input_axes, t, path) + Construct a Zarr dataset for the level n of a pyramid for the dimensions `input_axes`. It sets the type to `t` and saves it to `path/n` """ -function output_zarr(n, input_axes, t, path) - aggdims = agg_axis.(input_axes, 2^n) - s = length.(aggdims) +function output_zarr(n, input_axes, t, path, spatialdims;chunksizes=nothing) + spatial_axes = DD.dims(input_axes,spatialdims) + aggdims = agg_axis.(spatial_axes, 2^n) + spatial_axes_new = DD.setdims(input_axes,aggdims) + s = length.(spatial_axes_new) z = Zeros(t, s...) - yax = YAXArray(aggdims, z) - chunked = setchunks(yax , (1024, 1024)) + yax = YAXArray(spatial_axes_new, z) + if chunksizes === nothing + dchunksizes = map(d->DD.rebuild(d,1),input_axes) + dchunksizes = DD.setdims(dchunksizes,map(d->DD.rebuild(d,1024),spatial_axes)) + chunksizes = map(i->i.val,dchunksizes) + end + chunked = setchunks(yax , chunksizes) # This assumes that there is only the spatial dimensions to save ds = to_dataset(chunked, ) dssaved = savedataset(ds; path, skeleton=true, driver=:zarr) @@ -391,7 +414,7 @@ end Compute the data of the pyramids of a given data cube `ras`. This returns the data of the pyramids and the dimension values of the aggregated axes. """ -function getpyramids(reducefunc, ras;recursive=true, tilesize=256) +function getpyramids(reducefunc, ras;recursive=true, tilesize=256, spatial_dims=SpatialDim) input_axes = DD.dims(ras) n_level = compute_nlevels(ras, tilesize) if iszero(n_level) @@ -403,7 +426,9 @@ function getpyramids(reducefunc, ras;recursive=true, tilesize=256) outtype = Base.infer_return_type(reducefunc, (Matrix{eltype(ras)},)) #outtype = Missing <: eltype(ras) ? Union{Missing, outtypefunc} : outtypefunc outmin = output_arrays(pyramid_sizes, outtype) - fill_pyramids(ras,outmin,reducefunc,recursive; threaded=true) + ispatial_dims = DD.dimnum(DD.dims(ras),spatial_dims) + + fill_pyramids(ras,outmin,reducefunc,recursive, ispatial_dims; threaded=true) outmin, pyramid_axes end diff --git a/test/runtests.jl b/test/runtests.jl index 30cf63a..422b20c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,6 +47,7 @@ end @test length(axis.scene.plots) == 2 end +#= building an RGB pyramid doesn't work, need to think more about it. @testitem "Pyramid building RGB eltype" begin using PyramidScheme: PyramidScheme as PS using DimensionalData @@ -54,9 +55,10 @@ end data = rand(RGB, 2000,2000) dd = DimArray(data, (X(1:2000), Y(1:2000))) pyramid = PS.Pyramid(dd) - - + @test pyramid isa PS.Pyramid + @test eltype(pyramid) isa RGB end +=# @testitem "Helper functions" begin using PyramidScheme: PyramidScheme as PS @@ -140,6 +142,28 @@ end end end +@testitem "Building pyramid with additional dimension" begin + # The aim of this test is to check whether we can build a pyramid from a data cube with an extra dimension. + # We will only build the pyramids on the spatial dimensions and keep the other dimensions as is. + using YAXArrays + using Zarr + using PyramidScheme + using DimensionalData + s = (2048, 1024,3) + a = ones(s) + yax = YAXArray((X(1.:size(a,1)),Y(1.:size(a,2)), Z(1:3)), a) + path = tempname() *".zarr" + savecube(yax, path) + pyr = buildpyramids(path, resampling_method=sum) + pyrdisk = Pyramid(path) + for p in [pyr, pyrdisk] + @test p isa Pyramid + @test length(dims(p)) == 3 + @test size(p.levels[end]) == (256,128,3) + @test p.levels[1][1,1,1] == 4 + end + +end @testitem "Map of pyramids" begin using DimensionalData pyr1 = Pyramid(fill(1,X(1:128),Y(1:128)), tilesize=16)