Only code reorganization (src/ directory)

src/Plot/Plot.jl

@@ -0,0 +1,8 @@
+module Plot
+
+include("_plot.jl")
+
+export
+    wplotdots, wplotim
+
+end

src/mod/Plot.jl → src/Plot/_plot.jl

@@ -1,16 +1,17 @@
-module Plot
-export wplotdots, wplotim
-using ..Util, ..WT, ..Transforms
-using LinearAlgebra
+using ..Util
+using ..WT
+using ..Transforms
+
+using LinearAlgebra: norm
 
 # PLOTTING UTILITIES
 
 # return levels and detail coefficient centers on the interval [0,r) above (>=) threshold t
 # as tuple (d,l)
 function wplotdots(x::AbstractVector, t::Real=0, r::Real=1)
     isdyadic(x) || throw(ArgumentError("array must be of dyadic size"))
-	n = length(x)
-	c = wcount(x, t, level=0)
+    n = length(x)
+    c = wcount(x, t, level=0)
     d = Vector{Float64}(undef, c)
     l = Vector{Int}(undef, c)
     range = 0:1/n:1-eps()
@@ -22,15 +23,15 @@ function wplotdots(x::AbstractVector, t::Real=0, r::Real=1)
         for j = 0:J-1
             rind = 2^(J-1-j):2^(J-j):n
             for i in 1:dyadicdetailn(j)
-                if abs(x[dyadicdetailindex(j,i)]) >= t
+                if abs(x[dyadicdetailindex(j, i)]) >= t
                     d[k] = range[rind[i]]
                     l[k] = j
                     k += 1
                 end
             end
         end
     end
-    return (d,l)
+    return (d, l)
 end
 
 # return a matrix of detail coefficient values where row j+1 is level j
@@ -42,11 +43,11 @@ function wplotim(x::AbstractVector)
 
     @inbounds begin
         for j = 0:J-1
-	        dr = dyadicdetailrange(j)
-	        m = 2^(J-j)
-	        for i in eachindex(dr)
-		        A[j+1,1+(i-1)*m:i*m] .= x[dr[i]]
-	        end
+            dr = dyadicdetailrange(j)
+            m = 2^(J - j)
+            for i in eachindex(dr)
+                A[j+1, 1+(i-1)*m:i*m] .= x[dr[i]]
+            end
         end
     end
     return A
@@ -55,54 +56,50 @@ end
 # return an array of scaled detail coefficients and unscaled scaling coefficients
 # ready to be plotted as an image
 function wplotim(x::AbstractArray, L::Integer, wt::Union{DiscreteWavelet,Nothing}=nothing;
-                wabs::Bool=true, power::Real=0.7, pnorm::Real=1)
+    wabs::Bool=true, power::Real=0.7, pnorm::Real=1)
     isdyadic(x) || throw(ArgumentError("array must be of dyadic size"))
     dim = ndims(x)
     (dim == 2 || dim == 3) || throw(ArgumentError("dimension $(dim) not supported"))
-    n = size(x,1)
-    cn = size(x,3)  # color dimension
-    n == size(x,2) || throw(ArgumentError("array must be square"))
+    n = size(x, 1)
+    cn = size(x, 3)  # color dimension
+    n == size(x, 2) || throw(ArgumentError("array must be square"))
     (cn == 1 || cn == 3) || throw(ArgumentError("third dimension $(cn)  not supported"))
     J = ndyadicscales(n)
-    nsc = 2^(J-L)
+    nsc = 2^(J - L)
 
     # do wavelet transform
     if wt != nothing
-        if size(x,3)>1
+        if size(x, 3) > 1
             x = dwtc(x, wt, L)
         else
             x = dwt(x, wt, L)
         end
     end
 
     # scaling coefs
-    scs = x[1:nsc,1:nsc,:]
+    scs = x[1:nsc, 1:nsc, :]
     scale01!(scs)
 
     # detail coefs
     xts = copy(x)
-    wabs && (broadcast!(abs,xts,xts))
-    xts[1:nsc,1:nsc,:] .= 0
+    wabs && (broadcast!(abs, xts, xts))
+    xts[1:nsc, 1:nsc, :] .= 0
     scale01!(xts)
-    for j=1:n, i=1:n
-        @inbounds xts[i,j,:] .= norm(xts[i,j,:],pnorm).^(power)
+    for j = 1:n, i = 1:n
+        @inbounds xts[i, j, :] .= norm(xts[i, j, :], pnorm) .^ (power)
     end
 
     # merge and reshape final image
     scale01!(xts)
-    xts[1:nsc,1:nsc,:] = scs
+    xts[1:nsc, 1:nsc, :] = scs
     return xts
 end
 # scale elements of z to the interval [0,1]
 function scale01!(z)
     mi = minimum(z)
     ma = maximum(z)
     for i in eachindex(z)
-        @inbounds z[i] = (z[i] - mi)/(ma-mi)
+        @inbounds z[i] = (z[i] - mi) / (ma - mi)
     end
     return z
 end
-
-
-
-end

src/Threshold/Threshold.jl

@@ -0,0 +1,27 @@
+module Threshold
+
+using LinearAlgebra: rmul!, norm
+using Statistics: median!
+
+
+include("_threshold.jl")
+include("basis_functions.jl")
+include("denoising.jl")
+include("entropy.jl")
+
+
+export
+    # threshold
+    threshold!, threshold, HardTH, SoftTH, SemiSoftTH, SteinTH,
+    BiggestTH, PosTH, NegTH,
+
+    # denoising
+    DNFT, VisuShrink, denoise, noisest,
+
+    # basis functions
+    matchingpursuit, bestbasistree,
+
+    # entropy
+    coefentropy, Entropy, ShannonEntropy, LogEnergyEntropy
+
+end

src/Threshold/_threshold.jl

@@ -0,0 +1,129 @@
+using ..Util
+using ..WT
+using ..Transforms
+
+
+# THRESHOLD TYPES AND FUNCTIONS
+
+abstract type THType end
+struct HardTH <: THType end
+struct SoftTH <: THType end
+struct SemiSoftTH <: THType end
+struct SteinTH <: THType end
+struct BiggestTH <: THType end
+struct PosTH <: THType end
+struct NegTH <: THType end
+
+const DEFAULT_TH = HardTH()
+
+# biggest m-term approximation (best m-term approximation for orthogonal transforms)
+# result is m-sparse
+function threshold!(x::AbstractArray{<:Number}, TH::BiggestTH, m::Int)
+    @assert m >= 0
+    n = length(x)
+    m > n && (m = n)
+    ind = sortperm(x, alg=QuickSort, by=abs)
+    @inbounds begin
+        for i = 1:n-m
+            x[ind[i]] = 0
+        end
+    end
+    return x
+end
+
+# hard
+function threshold!(x::AbstractArray{<:Number}, TH::HardTH, t::Real)
+    @assert t >= 0
+    @inbounds begin
+        for i in eachindex(x)
+            if abs(x[i]) <= t
+                x[i] = 0
+            end
+        end
+    end
+    return x
+end
+
+# soft
+function threshold!(x::AbstractArray{<:Number}, TH::SoftTH, t::Real)
+    @assert t >= 0
+    @inbounds begin
+        for i in eachindex(x)
+            sh = abs(x[i]) - t
+            if sh < 0
+                x[i] = 0
+            else
+                x[i] = sign(x[i]) * sh
+            end
+        end
+    end
+    return x
+end
+
+# semisoft
+function threshold!(x::AbstractArray{<:Number}, TH::SemiSoftTH, t::Real)
+    @assert t >= 0
+    @inbounds begin
+        for i in eachindex(x)
+            if x[i] <= 2 * t
+                sh = abs(x[i]) - t
+                if sh < 0
+                    x[i] = 0
+                elseif sh - t < 0
+                    x[i] = sign(x[i]) * sh * 2
+                end
+            end
+        end
+    end
+    return x
+end
+
+# stein
+function threshold!(x::AbstractArray{<:Number}, TH::SteinTH, t::Real)
+    @assert t >= 0
+    @inbounds begin
+        for i in eachindex(x)
+            sh = 1 - t * t / (x[i] * x[i])
+            if sh < 0
+                x[i] = 0
+            else
+                x[i] = x[i] * sh
+            end
+        end
+    end
+    return x
+end
+
+# shrink negative elements to 0
+function threshold!(x::AbstractArray{<:Number}, TH::NegTH)
+    @inbounds begin
+        for i in eachindex(x)
+            if x[i] < 0
+                x[i] = 0
+            end
+        end
+    end
+    return x
+end
+
+# shrink positive elements to 0
+function threshold!(x::AbstractArray{<:Number}, TH::PosTH)
+    @inbounds begin
+        for i in eachindex(x)
+            if x[i] > 0
+                x[i] = 0
+            end
+        end
+    end
+    return x
+end
+
+# the non inplace functions
+function threshold(x::AbstractArray{T}, TH::THType, t::Real) where {T<:Number}
+    y = Vector{T}(undef, size(x))
+    return threshold!(copyto!(y, x), TH, t)
+end
+function threshold(x::AbstractArray{T}, TH::THType) where {T<:Number}
+    y = Vector{T}(undef, size(x))
+    return threshold!(copyto!(y, x), TH)
+end

src/Threshold/basis_functions.jl

@@ -0,0 +1,56 @@
+
+# BASIS FUNCTIONS
+
+# Matching Pursuit
+# see: Mallat (2009) p.642 "A wavelet tour of signal processing"
+# find sparse vector y such that ||x - f(y)|| < tol approximately
+# f is the operation of a M by N (M<N) dictionary/matrix
+# ft is a function defining the transpose of f
+function matchingpursuit(x::AbstractVector, f::Function, ft::Function, tol::Real, nmax::Int=-1, oop::Bool=false, N::Int=0)
+    @assert nmax >= -1
+    @assert tol > 0
+    r = x
+    n = 1
+
+    if !oop
+        y = zeros(eltype(x), length(ft(x)))
+    else # out of place functions f and ft
+        y = zeros(eltype(x), N)
+        tmp = similar(x, N)
+        ftr = similar(x, N)
+        aphi = similar(x, length(x))
+    end
+    spat = zeros(eltype(x), length(y))  # sparse for atom computation
+    nmax == -1 && (nmax = length(y))
+
+    while norm(r) > tol && n <= nmax
+        # find largest inner product
+        !oop && (ftr = ft(r))
+        oop && ft(ftr, r, tmp)
+        i = findmaxabs(ftr)
+
+        # project on i-th atom
+        spat[i] = ftr[i]
+        !oop && (aphi = f(spat))
+        oop && f(aphi, spat, tmp)
+        spat[i] = 0
+
+        # update residual, r = r - aphi
+        broadcast!(-, r, r, aphi)
+
+        y[i] += ftr[i]
+        n += 1
+    end
+    return y
+end
+function findmaxabs(x::AbstractVector)
+    m = abs(x[1])
+    k = 1
+    @inbounds for i in eachindex(x)
+        if abs(x[i]) > m
+            k = i
+            m = abs(x[i])
+        end
+    end
+    return k
+end