updates

Author: jason-yoon21

.gitignore

@@ -2,4 +2,5 @@
 test/.ipynb_checkpoints
 test/*.ipynb
 Manifest.toml
-Manifest-v1.11.toml
+Manifest-v1.11.toml
+*.txt

README.md

@@ -121,6 +121,31 @@ iphlct = iphlct2d(phlct, n)
 
 Reference: [Professor Saito's paper][paper] page 30 Chapter 6.2.2.
 
+## Examples
+
+Basic DST/IDST round-trip:
+
+```julia
+using .PolyHarmonicTrigTransforms
+
+x = randn(8)
+y = dst(x)
+x_rec = idst(y)
+@assert maximum(abs.(x_rec - 4 .* x ./ 4)) < 1e-12 # numeric check
+```
+
+PHLCT forward/backward round-trip (block size `N`):
+
+```julia
+using .PolyHarmonicTrigTransforms
+
+n = 8
+img = randn(n, n)
+coeffs = phlct_forward(img, n)
+img_rec = phlct_backward(coeffs, n)
+@assert norm(img - img_rec) / norm(img) < 1e-12
+```
+
 #### Image Testing
 
 To test images, first take your chosen image and convert it into a square image (i.e. size n x n where n is the new size of the image after using tools to square it).

src/dst.jl

@@ -30,23 +30,83 @@
 # See also: IDST
 
 module DST
-    using AbstractFFTs, FFTW
+    function _ensure_fftw()
+        if !isdefined(@__MODULE__, :FFTW)
+            Core.eval(Main, :(using FFTW))
+            Core.eval(@__MODULE__, :(const FFTW = Main.FFTW))
+        end
+        return nothing
+    end
 
-    export dst, dst_old
+    export dst, dst!, dst_u, plan_dst, dst_old
+
+    """
+    dst(x; dims=1)
+
+    Compute the Type-I Discrete Sine Transform (DST) of `x` along `dims`.
 
+    Arguments
+    - `x`: an `AbstractArray` (vector or matrix).
+    - `dims`: dimension to transform (default `1`).
+
+    Returns
+    - transformed array of the same shape as `x`.
+
+    Example
+    ```julia
+    y = dst([1.0,2.0,3.0])
+    ```
+    """
     function dst(x::AbstractArray, dims=1)
+        _ensure_fftw()
         N = size(x, dims)
         return FFTW.r2r(x, FFTW.RODFT00, dims) / 2
     end
 
     function dst_u(x::AbstractArray, dims=1)
+        _ensure_fftw()
         N = size(x, dims)
         new_x = FFTW.r2r(x, FFTW.RODFT00, dims)
 
         return new_x / (sqrt(2*N)) #unitory value
     end
 
+    """
+    plan_dst(x; dims=1, flags=FFTW.MEASURE)
+
+    Create an FFTW r2r plan for the DST (RODFT00) on `x`. The returned plan
+    can be executed with `plan * x`.
+    """
+    function plan_dst(x::AbstractArray; dims=1, flags=FFTW.MEASURE)
+        _ensure_fftw()
+        return FFTW.r2r_plan(x, FFTW.RODFT00, dims; flags=flags)
+    end
+
+    """
+    dst!(x; dims=1)
+
+    In-place DST that overwrites `x` with its transform when possible. Falls
+    back to a safe out-of-place compute and copy if the in-place routine is
+    not available on the current platform.
+    """
+    function dst!(x::AbstractArray; dims=1)
+        _ensure_fftw()
+        try
+            # Try to use FFTW in-place r2r if available
+            FFTW.r2r!(x, FFTW.RODFT00, dims)
+            x .*= 0.5
+            return x
+        catch _
+            # Fallback: compute out-of-place and copy into x
+            tmp = FFTW.r2r(x, FFTW.RODFT00, dims)
+            tmp .*= 0.5
+            x .= tmp
+            return x
+        end
+    end
+
     function dst_old(a::AbstractArray, n=nothing)
+        _ensure_fftw()
         if minimum(size(a)) == 1
             if size(a, 2) > 1
                 do_trans = 1
@@ -73,7 +133,7 @@ module DST
         y = zeros(2 * (n + 1), m)
         y[2:n+1, :] = aa
         y[n+3:2*(n+1), :] = -reverse(aa, dims=1)
-        yy = fft(y, (1,))
+        yy = FFTW.fft(y, (1,))
         b = yy[2:n+1, :] / (-2 * sqrt(-1 + 0im))
 
 

src/idst.jl

@@ -23,11 +23,31 @@
 # See also: DST
 # 
 module IDST
-
     include("dst.jl")
     using .DST
 
-    export idst, idst_old
+    # Lazily load FFTW when needed by IDST helpers.
+    function _ensure_fftw()
+        if !isdefined(@__MODULE__, :FFTW)
+            Core.eval(Main, :(using FFTW))
+            Core.eval(@__MODULE__, :(const FFTW = Main.FFTW))
+        end
+        return nothing
+    end
+
+    export idst, idst!, plan_idst, idst_old
+    """
+    idst(y; dims=1)
+
+    Inverse Type-I Discrete Sine Transform (IDST) of `y` along `dims`.
+
+    Arguments
+    - `y`: an `AbstractArray` of DST coefficients.
+    - `dims`: dimension to transform (default `1`).
+
+    Returns
+    - reconstructed array in the original domain.
+    """
     function idst(y::AbstractArray, dims=1)
 
         N = size(y, dims)
@@ -37,6 +57,36 @@ module IDST
         return x * 4
     end
 
+    """
+    plan_idst(y; dims=1, flags=FFTW.MEASURE)
+
+    Create an FFTW r2r plan suitable for the IDST (uses the DST plan internally).
+    """
+    function plan_idst(y::AbstractArray; dims=1, flags=FFTW.MEASURE)
+        _ensure_fftw()
+        return FFTW.r2r_plan(y, FFTW.RODFT00, dims; flags=flags)
+    end
+
+    """
+    idst!(y; dims=1)
+
+    In-place inverse DST: overwrite `y` with the reconstructed data when
+    possible. Falls back to computing out-of-place and copying the result.
+    """
+    function idst!(y::AbstractArray; dims=1)
+        _ensure_fftw()
+        try
+            # In-place inverse DST: perform r2r in-place then scale by 1/(N+1)
+            FFTW.r2r!(y, FFTW.RODFT00, dims)
+            y ./= (size(y, dims) + 1)
+            return y
+        catch _
+            tmp = idst(y, dims)
+            y .= tmp
+            return y
+        end
+    end
+
     function idst_old(y::AbstractArray, n=nothing, dims=1)
         if n === nothing
             n = size(y, 1)

src/llst.jl

@@ -23,6 +23,7 @@ module LLST
     using .SOLVELAPLACE
 
     export llst2d, illst2d
+    export llst2d!, illst2d!
 
     # LLFT2D 2D Laplace Local Fourier Transform
     # 
@@ -51,11 +52,24 @@ module LLST
         inIX1 = 2:m-1
         inIX2 = 2:n-1
         rnorm = sqrt(4 / (m - 1) / (n - 1))
-        interior = data[inIX1, inIX2] .- solvelaplace(data)
+        sl = solvelaplace(data)
+        interior = data[inIX1, inIX2] .- sl[inIX1, inIX2]
         trans[inIX1, inIX2] = rnorm .* dst(dst(interior)')'
         return trans
     end
 
+    """
+    llst2d!(data)
+
+    In-place variant of `llst2d` that overwrites `data` with the transform
+    result. Falls back to out-of-place computation and broadcasting assignment.
+    """
+    function llst2d!(data::AbstractVecOrMat)
+        res = llst2d(data)
+        data .= res
+        return data
+    end
+
     # ILLFT2D 2D Inverse Laplace Local Fourier Transform
     # 
     # Applies the inverse LLFT transform to image (no blocking). Note that 
@@ -94,10 +108,23 @@ module LLST
 
         data = float(data)
         data[inIX1[:], inIX2[:]] = rnorm .* dst(dst(data[inIX1[:], inIX2[:]])')'
-        trans[inIX1[:], inIX2[:]] = data[inIX1[:], inIX2[:]] .+ solvelaplace(trans)
+        sl = solvelaplace(trans)
+        trans[inIX1[:], inIX2[:]] = data[inIX1[:], inIX2[:]] .+ sl[inIX1[:], inIX2[:]]
 
         return trans
     end
 
+    """
+    illst2d!(data)
+
+    In-place variant of `illst2d` that overwrites `data` with the inverse
+    transform result.
+    """
+    function illst2d!(data::AbstractVecOrMat)
+        res = illst2d(data)
+        data .= res
+        return data
+    end
+
 
 end

src/llst2.jl

@@ -10,6 +10,7 @@ module LLST2
     using .SOLVELAPLACE
 
     export llst2
+    export llst2!
 
     function llst2(data::AbstractVecOrMat, ll::AbstractArray, inverse::Bool=false)
         if (inverse != false)
@@ -34,6 +35,18 @@ module LLST2
         end
     end
 
+    """
+    llst2!(data, ll; inverse=false)
+
+    In-place variant of `llst2` that overwrites `data` with the transformed
+    result. Uses broadcasting assignment as a safe fallback.
+    """
+    function llst2!(data::AbstractVecOrMat, ll::AbstractArray, inverse::Bool=false)
+        res = llst2(data, ll, inverse)
+        data .= res
+        return data
+    end
+
 
     # Function for forward transform of boundary
     function fbdrytrans(bdry::AbstractVecOrMat)