Question : Working on rectangular matrix -> multilevel fully separable decomposition

[aTrotier]

I am working with non-squared matrix 2D (x != y) or 3D (x != y != z). For example an image of 220 x 160, when applying the DWT it apply a DWT of level 2 (restricted by the size 220)

1. Is that possible to perform a multilevel fully separable decomposition WT like : https://pywavelets.readthedocs.io/en/latest/ref/nd-dwt-and-idwt.html#multilevel-fully-separable-decomposition-fswavedecn I did not find the function for that but maybe I am missing something :)

2. An other approach will be to pad the data. In that case does the operator is still orthogonal (Adjoint = Inverse)

Thanks :)

[gummif]

I don´t think there is. You could do it yourself if I understand correctly, just transform fully along one dimension then again fully the other direction and keep track of everything yourself.

[aTrotier]

I have a working example for that :

using ImageUtils
using Wavelets
using Plots

ph = shepp_logan(128)
heatmap(ph)

W_ph = dwt(ph,wavelet(WT.db2))
heatmap(abs.(W_ph),clim=(0,1))


## perform Wavelet on rectangular matrix with different maximum deph level

ph2 = ph[:,1:80]
sph2 = size(ph2)
MaxL = Int[]
for (i,val) in enumerate(sph2)
    push!(MaxL,maxtransformlevels(val))
end
display(MaxL)

# when using dwt only max level transform is used -> 4
W_ph = dwt(ph2,wavelet(WT.db2))
heatmap(abs.(W_ph),clim=(0,1))

# we can perform dwt along one axis for max transform = 7
res = zeros(Float64,sph2)
for i in 1:sph2[1]
    res[i,:] = dwt(ph2[i,:],wavelet(WT.db2))
end

# and then the second axis for maxtransform = 4
for j in 1:sph2[2]
    res[:,j] = dwt(res[:,j],wavelet(WT.db2))
end
heatmap(res)

## get back the original image
im = zeros(Float64,sph2)
for j in 1:sph2[2]
    im[:,j] = idwt(res[:,j],wavelet(WT.db2))
end

for i in 1:sph2[1]
    im[i,:] = idwt(im[i,:],wavelet(WT.db2))
end
heatmap(im)

Do you want to implement a function for nd-dwt/idwt like nd_dwt(x::Array{T,D},wt::OrthoFilter [;L ,dims]) where dwt is fully transformed only along dims ?

I can implement that only in MRIReco but maybe @JeffFessler will also be interested to use that implementation

---

Related to https://github.com/tknopp/SparsityOperators.jl/issues/12

[JeffFessler]

There could be many cases with non-square multi-dimensional data where users could want different levels across different dimensions. So I agree that it would be desirable to have that functionality here in Wavelets.jl.

One could generalize dwt and idwt to allow the L argument to be a NTuple{D,Int} where x::AbstractArray{T,D}. The catch is that using, say, L=(4,4) would not produce the same results as L=4 because the former would use a separable decomposition whereas the latter would not (I think). I wonder if there is a more unified way.

[aTrotier]

Oh, I thought that results will be the same but your are right, it is totally different :

Do you think it is better for CS-MRI to use the standard dwt ?

@uecker what dwt transform are you using in BART ?

[uecker]

I guess I implemented a more unified way: The hierarchical decomposition works correctly for arbitrary dimensions with different number of levels by recursing into the low-pass segment.

[aTrotier]

@uecker something like that :

## try a more unified way
ph = shepp_logan(128)
ph2 = ph[:,1:100]
sph2 = size(ph2)
MaxL = Int[]
for (i,val) in enumerate(sph2)
    push!(MaxL,maxtransformlevels(val))
end
display(MaxL)
L = minimum(MaxL)
W_tmp = dwt(ph2,wavelet(WT.db2),L)
p1 = heatmap(abs.(W_tmp),title ="L = $L - standard",clim=(0,1),aspect_ratio = 1,legend = :none , axis=nothing)

# remove
W_lp = W_tmp[1:Int(sph2[1]/2^L),1:Int(sph2[2]/2^L)]
heatmap(abs.(W_lp))

tmp = similar(W_lp)
@views for j in 1:size(W_lp)[2]
    dwt!(tmp[:,j], W_lp[:,j],wavelet(WT.db2))
 end

 heatmap(abs.(tmp))

 W_tmp[1:Int(sph2[1]/2^L),1:Int(sph2[2]/2^L)] .= tmp
 p2 = heatmap(abs.(W_tmp),title ="L = $L + supp dwt along dims",clim=(0,1),aspect_ratio = 1,legend = :none , axis=nothing)

 plot(p1,p2)

[JeffFessler]

That looks better to me. BTW to simplify the code I think you could use broadcast: MaxL = maxtransformlevels.(size(ph2))

Another comment is that I never threshold the approximation coefficients (low-pass segment) because it is not expected to be sparse. I am not sure how well known or used that variation is in MRI. Here is an example: https://juliaimagerecon.github.io/Examples/generated/mri/2-cs-wl-l1-2d/#Wavelet-sparsity-in-synthesis-form I guess that remark is more relevant to MRreco than here...

[tknopp]

I guess that remark is more relevant to MRreco than here...

I don't think so. Thresholding is part of Wavelets.jl (https://github.com/JuliaDSP/Wavelets.jl#thresholding). So it would indeed be interesting what is implemented there.