Open aTrotier opened 2 years ago
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.
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
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.
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 ?
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.
@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)
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...
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.
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 :)