From 3c50c6ce77bed1dfd432c8be2ac3f76919d483a9 Mon Sep 17 00:00:00 2001 From: miladbader Date: Tue, 17 Mar 2026 15:08:56 +0000 Subject: [PATCH 1/4] generalize JopSlantStack3D to handle 5D input (with xy passive dimensions) --- Project.toml | 2 +- src/JetPackTransforms.jl | 2 +- src/jop_slantstack.jl | 221 ++++++++++++++++++++++++++++++--------- test/jop_slantstack.jl | 31 ++++++ 4 files changed, 202 insertions(+), 54 deletions(-) diff --git a/Project.toml b/Project.toml index 8d45e28..96a1588 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "JetPackTransforms" uuid = "e0630bc5-6e1f-509c-a3e4-0aea24bd2b1c" -version = "0.3.0" +version = "0.4.0" [deps] DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" diff --git a/src/JetPackTransforms.jl b/src/JetPackTransforms.jl index 64dc618..66393ac 100644 --- a/src/JetPackTransforms.jl +++ b/src/JetPackTransforms.jl @@ -6,6 +6,6 @@ include("jop_dct.jl") include("jop_dwt.jl") include("jop_fft.jl") include("jop_sft.jl") -include("jop_slantstack.jl") # requires JopTaper from JetPack and WaveFD.shiftforward! and WaveFD.shiftadjoint! from WaveFD.jl +include("jop_slantstack.jl") # requires JopTaper from JetPack end diff --git a/src/jop_slantstack.jl b/src/jop_slantstack.jl index d1ba160..ff9a8c7 100644 --- a/src/jop_slantstack.jl +++ b/src/jop_slantstack.jl @@ -171,9 +171,10 @@ end A = JopSlantStack3D(dom[; dz=1.0, dhx=1.0, dhy=1.0, hx0=0.0, hy0=0.0, ...]) where `A` is the 3D slant-stack operator mapping for `z-hx-hy` to `z-θ-ϕ` (depth mode) or `t-hx-hy` to `tau-px-py` (time mode). -The domain of the operator is `nz` x `nhx` x `nhy` with precision T, `dz` is the depth spacing (or time interval), -`dhx` and `dhy` are the offset spacings, and `hx0` and `hy0` are the origins of the offset axes. The additional named optional arguments -along with their default values are, +The domain of the operator is typically `nz` x `nhx` x `nhy` with precision T, `dz` is the depth spacing (or time interval), +`dhx` and `dhy` are the offset spacings, and `hx0` and `hy0` are the origins of the offset axes. +The domain can also be `nz` x `ny` x `nx` x `nhx` x `nhy` where the `y` and `x` dimensions are passive. +The additional named optional arguments along with their default values are, * `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-hx-hy` or `t-hx-hy`. * `theta=-45:1.0:45` - range of incidence angles used when `mode="depth"`. @@ -227,7 +228,16 @@ function JopSlantStack3D( dhy < 0.0 && error("expected dhy>0.0, got dhy=$(dhy)") mode == "depth" && any(abs.(theta) .>= 90.0) && error("for mode='depth', all theta values must be < 90 degrees in absolute value") - nz,nhx,nhy = size(dom) + ndim = length(size(dom)) + if ndim == 3 + nz,nhx,nhy = size(dom) + active_dims = (1,2,3) + elseif ndim == 5 + nz,ny,nx,nhx,nhy = size(dom) + active_dims = (1,4,5) + else + error("unsupported domain size: $(size(dom))") + end # kz nzfft = nextprod([2,3,5,7], round(Int, nz*(1 + padz))) @@ -266,10 +276,10 @@ function JopSlantStack3D( azimuth = T(azimuth) # tapers - TX = JopTaper(dom, (1,2,3), (taperz[1],taperhx[1],taperhy[1]), (taperz[2],taperhx[2],taperhy[2])) - TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1, nhxfft, nhyfft), (1,2,3), (taperkz[1], taperkhx[1], taperkhy[1]), (taperkz[2], taperkhx[2], taperkhy[2]), mode=(:normal,:fftshift,:fftshift)) + TX = JopTaper(dom, active_dims, (taperz[1],taperhx[1],taperhy[1]), (taperz[2],taperhx[2],taperhy[2])) + TK = JopTaper(ndim == 3 ? JetSpace(Complex{eltype(dom)}, div(nzfft,2)+1, nhxfft, nhyfft) : JetSpace(Complex{eltype(dom)}, div(nzfft,2)+1, ny, nx, nhxfft, nhyfft), active_dims, (taperkz[1], taperkhx[1], taperkhy[1]), (taperkz[2], taperkhx[2], taperkhy[2]), mode=(:normal,:fftshift,:fftshift)) - JopLn(dom = dom, rng = JetSpace(T, nz, length(px), length(py)), df! = JopSlantStack3D_df!, df′! = JopSlantStack3D_df′!, + JopLn(dom = dom, rng = ndim == 3 ? JetSpace(T, nz, length(px), length(py)) : JetSpace(T, nz, ny, nx, length(px), length(py)), df! = JopSlantStack3D_df!, df′! = JopSlantStack3D_df′!, s = (;mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK)) end export JopSlantStack3D @@ -293,37 +303,89 @@ function JopSlantStack3D_df!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; mode, is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft) D = zeros(eltype(M), size(M,1), npx, npy) - @threads for ipx = 1:npx - for ipy = 1:npy - num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) - denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) - for ikz = 1:div(nzfft,2)+1 - ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) + @threads for ip = 1:npx*npy + ipx = div(ip-1, npy) + 1 + ipy = mod(ip-1, npy) + 1 + num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) + denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) + for ikz = 1:div(nzfft,2)+1 + ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) + + is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue + + if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) + D[ikz,ipx,ipy] = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + continue + end - is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue + ax_m1 = (khx[ikhx_p1] - _khx)/dhx + ax_p1 = 1 - ax_m1 + ay_m1 = (khy[ikhy_p1] - _khy)/dhy + ay_p1 = 1 - ay_m1 - if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) - D[ikz,ipx,ipy] = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) - continue - end + d_p1_p1 = M[ikz,ikhx_p1,ikhy_p1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + d_p1_m1 = M[ikz,ikhx_p1,ikhy_m1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) + d_m1_p1 = M[ikz,ikhx_m1,ikhy_p1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) + d_m1_m1 = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + + D[ikz,ipx,ipy] = ax_m1*ay_m1*d_m1_m1 + ax_m1*ay_p1*d_m1_p1 + ax_p1*ay_m1*d_p1_m1 + ax_p1*ay_p1*d_p1_p1 + end + end + + _d = brfft(D, nzfft, 1) + d .= _d[1:nz,1:npx,1:npy] ./ nzfft +end + +function JopSlantStack3D_df!(d::AbstractArray{T,5}, m::AbstractArray{T,5}; mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK, kwargs...) where {T} + ny, nx, nhx, nhy, npx, npy, dhx, dhy = size(m,2), size(m,3), size(m,4), size(m,5), length(px), length(py), abs(khx[2]-khx[1]), abs(khy[2]-khy[1]) + nz = size(d, 1) + + mpad = zeros(T, nzfft, ny, nx, nhxfft, nhyfft) + mpad[1:nz,1:ny,1:nx,1:nhx,1:nhy] = TX*m + M = TK*rfft(mpad, (1, 4, 5)) + + if mode == "time" + compute_kh = slantstack_khxy_from_frequency + elseif dip == 0 + compute_kh = slantstack_khxy_from_kz + else + compute_kh = slantstack_khxy_from_kz_geologic + end - ax_m1 = (khx[ikhx_p1] - _khx)/dhx - ax_p1 = 1 - ax_m1 - ay_m1 = (khy[ikhy_p1] - _khy)/dhy - ay_p1 = 1 - ay_m1 + is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft) + + D = zeros(eltype(M), size(M,1), ny, nx, npx, npy) + @threads for ip = 1:npx*npy + ipx = div(ip-1, npy) + 1 + ipy = mod(ip-1, npy) + 1 + num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) + denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) + for ikz = 1:div(nzfft,2)+1 + ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) - d_p1_p1 = M[ikz,ikhx_p1,ikhy_p1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) - d_p1_m1 = M[ikz,ikhx_p1,ikhy_m1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) - d_m1_p1 = M[ikz,ikhx_m1,ikhy_p1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) - d_m1_m1 = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue - D[ikz,ipx,ipy] = ax_m1*ay_m1*d_m1_m1 + ax_m1*ay_p1*d_m1_p1 + ax_p1*ay_m1*d_p1_m1 + ax_p1*ay_p1*d_p1_p1 + if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) + @views D[ikz,:,:,ipx,ipy] .= M[ikz,:,:,ikhx_m1,ikhy_m1] .* exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + continue end + + ax_m1 = (khx[ikhx_p1] - _khx)/dhx + ax_p1 = 1 - ax_m1 + ay_m1 = (khy[ikhy_p1] - _khy)/dhy + ay_p1 = 1 - ay_m1 + + d_p1_p1 = M[ikz,:,:,ikhx_p1,ikhy_p1] .* exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + d_p1_m1 = M[ikz,:,:,ikhx_p1,ikhy_m1] .* exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) + d_m1_p1 = M[ikz,:,:,ikhx_m1,ikhy_p1] .* exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) + d_m1_m1 = M[ikz,:,:,ikhx_m1,ikhy_m1] .* exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + + @views D[ikz,:,:,ipx,ipy] .= (ax_m1*ay_m1) .* d_m1_m1 .+ (ax_m1*ay_p1) .* d_m1_p1 .+ (ax_p1*ay_m1) .* d_p1_m1 .+ (ax_p1*ay_p1) .* d_p1_p1 end end _d = brfft(D, nzfft, 1) - d .= _d[1:nz,1:npx,1:npy] ./ nzfft + d .= _d[1:nz,:,:,1:npx,1:npy] ./ nzfft end function JopSlantStack3D_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK, kwargs...) where {T} @@ -346,39 +408,94 @@ function JopSlantStack3D_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; mo is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft) M = zeros(Complex{T}, div(nzfft,2)+1, nhxfft, nhyfft) - for ipx = 1:npx - for ipy = 1:npy - num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) - denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) - for ikz = 1:div(nzfft,2)+1 - ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) + for ip = 1:npx*npy + ipx = div(ip-1, npy) + 1 + ipy = mod(ip-1, npy) + 1 + num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) + denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) + for ikz = 1:div(nzfft,2)+1 + ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) + + is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue + + if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) + M[ikz,ikhx_p1,ikhy_p1] += D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + continue + end - is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue + ax_m1 = (khx[ikhx_p1] - _khx)/dhx + ax_p1 = 1 - ax_m1 + ay_m1 = (khy[ikhy_p1] - _khy)/dhy + ay_p1 = 1 - ay_m1 - if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) - M[ikz,ikhx_p1,ikhy_p1] += D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) - continue - end + m_p1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + m_p1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) + m_m1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) + m_m1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + + M[ikz,ikhx_p1,ikhy_p1] += ax_p1*ay_p1*m_p1_p1 + M[ikz,ikhx_p1,ikhy_m1] += ax_p1*ay_m1*m_p1_m1 + M[ikz,ikhx_m1,ikhy_p1] += ax_m1*ay_p1*m_m1_p1 + M[ikz,ikhx_m1,ikhy_m1] += ax_m1*ay_m1*m_m1_m1 + end + end + + m .= TX * (brfft(TK * M, nzfft)[1:nz,1:nhx,1:nhy]) +end + +function JopSlantStack3D_df′!(m::AbstractArray{T,5}, d::AbstractArray{T,5}; mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK, kwargs...) where {T} + ny, nx, nhx, nhy, npx, npy, dhx, dhy = size(m,2), size(m,3), size(m,4), size(m,5), length(px), length(py), abs(khx[2]-khx[1]), abs(khy[2]-khy[1]) + nz = size(d, 1) + + dpad = zeros(T, nzfft, ny, nx, npx, npy) + dpad[1:nz,:,:,:,:] = d - ax_m1 = (khx[ikhx_p1] - _khx)/dhx - ax_p1 = 1 - ax_m1 - ay_m1 = (khy[ikhy_p1] - _khy)/dhy - ay_p1 = 1 - ay_m1 + D = (rfft(dpad, 1) ./ nzfft) + + if mode == "time" + compute_kh = slantstack_khxy_from_frequency + elseif dip == 0 + compute_kh = slantstack_khxy_from_kz + else + compute_kh = slantstack_khxy_from_kz_geologic + end + + is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft) + + M = zeros(Complex{T}, div(nzfft,2)+1, ny, nx, nhxfft, nhyfft) + for ip = 1:npx*npy + ipx = div(ip-1, npy) + 1 + ipy = mod(ip-1, npy) + 1 + num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy]) + denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2) + for ikz = 1:div(nzfft,2)+1 + ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) - m_p1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) - m_p1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) - m_m1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) - m_m1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue - M[ikz,ikhx_p1,ikhy_p1] += ax_p1*ay_p1*m_p1_p1 - M[ikz,ikhx_p1,ikhy_m1] += ax_p1*ay_m1*m_p1_m1 - M[ikz,ikhx_m1,ikhy_p1] += ax_m1*ay_p1*m_m1_p1 - M[ikz,ikhx_m1,ikhy_m1] += ax_m1*ay_m1*m_m1_m1 + if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1) + @views M[ikz,:,:,ikhx_p1,ikhy_p1] .+= D[ikz,:,:,ipx,ipy] .* exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + continue end + + ax_m1 = (khx[ikhx_p1] - _khx)/dhx + ax_p1 = 1 - ax_m1 + ay_m1 = (khy[ikhy_p1] - _khy)/dhy + ay_p1 = 1 - ay_m1 + + m_p1_p1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) + m_p1_m1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) + m_m1_p1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) + m_m1_m1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) + + M[ikz,:,:,ikhx_p1,ikhy_p1] .+= (ax_p1*ay_p1) .* m_p1_p1 + M[ikz,:,:,ikhx_p1,ikhy_m1] .+= (ax_p1*ay_m1) .* m_p1_m1 + M[ikz,:,:,ikhx_m1,ikhy_p1] .+= (ax_m1*ay_p1) .* m_m1_p1 + M[ikz,:,:,ikhx_m1,ikhy_m1] .+= (ax_m1*ay_m1) .* m_m1_m1 end end - m .= TX * (brfft(TK * M, nzfft)[1:nz,1:nhx,1:nhy]) + m .= TX * (brfft(TK * M, nzfft, (1,4,5))[1:nz,:,:,1:nhx,1:nhy]) end @inline function slantstack_khxy_from_kz(ikz::Int64, ipx::Int64, ipy::Int64, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) diff --git a/test/jop_slantstack.jl b/test/jop_slantstack.jl index e88ae39..9270abe 100644 --- a/test/jop_slantstack.jl +++ b/test/jop_slantstack.jl @@ -118,6 +118,11 @@ end lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A))) @test isapprox(lhs,rhs,rtol=1e-4) + + A = JopSlantStack3D(JetSpace(T, 32, 3, 4, 8, 5); dz=10.0, dhx=10.0, dhy=5.0, hx0=-1000.0, hy0=-500.0, mode=mode, dip=dip, azimuth=azimuth, taperz=(0.3,0.3), taperhx=(0.3,0.3), taperhy=(0.3,0.3), taperkz=(0.3,0.3), taperkhx=(0.3,0.3), taperkhy=(0.3,0.3)) + + lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A))) + @test isapprox(lhs,rhs,rtol=1e-4) end @testset "JetSlantStack3D, time-depth parity" begin @@ -150,6 +155,32 @@ end @test err_d2 < 1e-7 end +@testset "JetSlantStack3D, 3D vs 5D parity" begin + theta = collect(-45:1.0:45) + phi = [0.0] + px = @. -tan(deg2rad(theta)) + py = [0.0] + + A1 = JopSlantStack3D(JetSpace(Float64, 64, 8, 8); mode="depth", theta=theta, phi = phi, px=px, py=py, dz=0.01, dhx=0.01, dhy=0.01, hx0=-0.64, hy0=-0.64) + A2 = JopSlantStack3D(JetSpace(Float64, 64, 1, 1, 8, 8); mode="depth", theta=theta, phi = phi, px=px, py=py, dz=0.01, dhx=0.01, dhy=0.01, hx0=-0.64, hy0=-0.64) + m = zeros(domain(A1)) + for i = 1:8 + m[div(i,4)+20,i,:] .= 1.0 + m[div(i,2)+20,:,i] .= 0.5 + end + m[16,4,4] = 1 + m[32,5,7] = 1 + + d1 = A1*m + m = reshape(m, size(domain(A2))) + d2 = A2*m + d2 = reshape(d2, size(d1)) + + err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1)) + @show err_d + @test err_d < 1e-7 +end + @testset "JetSlantStackShiftSum3D, dot product test, T = $(T), mode = $(mode), dip = $(dip)" for T in (Float32, Float64), mode in ("depth", "time"), (dip, azimuth) in ((0.0,0.0), (30.0, 60.0), ([0.0], [0.0])) if isa(dip, Array) dip = 90 .* rand(32) From bb51faf619579e754e72d3f3ff1a789852a2c756 Mon Sep 17 00:00:00 2001 From: miladbader Date: Tue, 17 Mar 2026 15:58:09 +0000 Subject: [PATCH 2/4] Add random seed for UT --- test/jop_slantstack.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/jop_slantstack.jl b/test/jop_slantstack.jl index 9270abe..7d18c37 100644 --- a/test/jop_slantstack.jl +++ b/test/jop_slantstack.jl @@ -1,4 +1,4 @@ -using JetPackTransforms, Jets, Test, LinearAlgebra +using JetPackTransforms, Jets, Test, LinearAlgebra, Random # depth mode @testset "JetSlantStack, dot product test" for T in (Float32, Float64) @@ -114,6 +114,7 @@ end end @testset "JetSlantStack3D, dot product test, T = $(T), mode = $(mode), dip = $(dip)" for T in (Float32, Float64), mode in ("depth", "time"), (dip, azimuth) in ((0.0,0.0), (30.0, 60.0)) + Random.seed!(1234) A = JopSlantStack3D(JetSpace(T, 32, 8, 5); dz=10.0, dhx=10.0, dhy=5.0, hx0=-1000.0, hy0=-500.0, mode=mode, dip=dip, azimuth=azimuth, taperz=(0.3,0.3), taperhx=(0.3,0.3), taperhy=(0.3,0.3), taperkz=(0.3,0.3), taperkhx=(0.3,0.3), taperkhy=(0.3,0.3)) lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A))) From cce98da250751067698cc1ed1e7ccc8573e468b4 Mon Sep 17 00:00:00 2001 From: miladbader Date: Tue, 17 Mar 2026 17:27:05 +0000 Subject: [PATCH 3/4] avoid array allocation in multithreading --- src/jop_slantstack.jl | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/jop_slantstack.jl b/src/jop_slantstack.jl index ff9a8c7..335927f 100644 --- a/src/jop_slantstack.jl +++ b/src/jop_slantstack.jl @@ -375,12 +375,10 @@ function JopSlantStack3D_df!(d::AbstractArray{T,5}, m::AbstractArray{T,5}; mode, ay_m1 = (khy[ikhy_p1] - _khy)/dhy ay_p1 = 1 - ay_m1 - d_p1_p1 = M[ikz,:,:,ikhx_p1,ikhy_p1] .* exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) - d_p1_m1 = M[ikz,:,:,ikhx_p1,ikhy_m1] .* exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) - d_m1_p1 = M[ikz,:,:,ikhx_m1,ikhy_p1] .* exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) - d_m1_m1 = M[ikz,:,:,ikhx_m1,ikhy_m1] .* exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) - - @views D[ikz,:,:,ipx,ipy] .= (ax_m1*ay_m1) .* d_m1_m1 .+ (ax_m1*ay_p1) .* d_m1_p1 .+ (ax_p1*ay_m1) .* d_p1_m1 .+ (ax_p1*ay_p1) .* d_p1_p1 + @views D[ikz,:,:,ipx,ipy] .= (ax_m1*ay_m1*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0))) .* M[ikz,:,:,ikhx_m1,ikhy_m1] .+ + (ax_m1*ay_p1*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0))) .* M[ikz,:,:,ikhx_m1,ikhy_p1] .+ + (ax_p1*ay_m1*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0))) .* M[ikz,:,:,ikhx_p1,ikhy_m1] .+ + (ax_p1*ay_p1*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0))) .* M[ikz,:,:,ikhx_p1,ikhy_p1] end end @@ -483,15 +481,10 @@ function JopSlantStack3D_df′!(m::AbstractArray{T,5}, d::AbstractArray{T,5}; mo ay_m1 = (khy[ikhy_p1] - _khy)/dhy ay_p1 = 1 - ay_m1 - m_p1_p1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0)) - m_p1_m1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0)) - m_m1_p1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0)) - m_m1_m1 = D[ikz,:,:,ipx,ipy].* exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0)) - - M[ikz,:,:,ikhx_p1,ikhy_p1] .+= (ax_p1*ay_p1) .* m_p1_p1 - M[ikz,:,:,ikhx_p1,ikhy_m1] .+= (ax_p1*ay_m1) .* m_p1_m1 - M[ikz,:,:,ikhx_m1,ikhy_p1] .+= (ax_m1*ay_p1) .* m_m1_p1 - M[ikz,:,:,ikhx_m1,ikhy_m1] .+= (ax_m1*ay_m1) .* m_m1_m1 + @views M[ikz,:,:,ikhx_p1,ikhy_p1] .+= (ax_p1*ay_p1*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0))) .* D[ikz,:,:,ipx,ipy] + @views M[ikz,:,:,ikhx_p1,ikhy_m1] .+= (ax_p1*ay_m1*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0))) .* D[ikz,:,:,ipx,ipy] + @views M[ikz,:,:,ikhx_m1,ikhy_p1] .+= (ax_m1*ay_p1*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0))) .* D[ikz,:,:,ipx,ipy] + @views M[ikz,:,:,ikhx_m1,ikhy_m1] .+= (ax_m1*ay_m1*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0))) .* D[ikz,:,:,ipx,ipy] end end From 14bcb5390645020fd00670fefad4a3273eab05ea Mon Sep 17 00:00:00 2001 From: miladbader Date: Wed, 18 Mar 2026 20:37:59 +0000 Subject: [PATCH 4/4] fix bug in geologic mapping function --- src/jop_slantstack.jl | 4 ++-- test/jop_slantstack.jl | 46 +++++++++++++++++++++++++++++++++++++++++- 2 files changed, 47 insertions(+), 3 deletions(-) diff --git a/src/jop_slantstack.jl b/src/jop_slantstack.jl index 335927f..222d9bf 100644 --- a/src/jop_slantstack.jl +++ b/src/jop_slantstack.jl @@ -509,8 +509,8 @@ end end @inline function slantstack_khxy_from_kz_geologic(ikz::Int64, ipx::Int64, ipy::Int64, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom) - _khx = kz[ikz]*px[ipx]*(cos(py[ipy])/denom + sin(py[ipy]*num/denom)) - _khy = kz[ikz]*px[ipx]*(sin(py[ipy])/denom - cos(py[ipy]*num/denom)) + _khx = kz[ikz]*px[ipx]*(cos(py[ipy]) + sin(py[ipy])*num)/denom + _khy = kz[ikz]*px[ipx]*(sin(py[ipy]) - cos(py[ipy])*num)/denom ikhx_m1 = floor(Int64, _khx/khx[2]) + 1 ikhx_p1 = ceil(Int64, _khx/khx[2]) + 1 diff --git a/test/jop_slantstack.jl b/test/jop_slantstack.jl index 7d18c37..09327b8 100644 --- a/test/jop_slantstack.jl +++ b/test/jop_slantstack.jl @@ -182,6 +182,28 @@ end @test err_d < 1e-7 end +@testset "JetSlantStack3D, dip invariance" begin + theta = collect(-45:1.0:45) + phi = [0.0] + + A1 = JopSlantStack3D(JetSpace(Float64, 128, 129, 5); mode="depth", theta=theta, phi=phi, dz=0.01, hx0=-0.64, dhx=0.01, hy0=0, dhy=0.01, dip=0.0, azimuth=0.0) + A2 = JopSlantStack3D(JetSpace(Float64, 128, 129, 5); mode="depth", theta=theta, phi=phi, dz=0.01, hx0=-0.64, dhx=0.01, hy0=0, dhy=0.01, dip=30.0, azimuth=0.0) + m = zeros(domain(A1)) + for i = 1:129 + m[div(i,4)+20,i,1] = 1.0 + m[div(i,3)+10,i,2] = -0.5 + end + m[16,64,1] = 1 + m[64,96,5] = -0.2 + d1 = A1*m + d2 = A2*m + + err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1)) + + @show err_d + @test err_d < 1e-7 +end + @testset "JetSlantStackShiftSum3D, dot product test, T = $(T), mode = $(mode), dip = $(dip)" for T in (Float32, Float64), mode in ("depth", "time"), (dip, azimuth) in ((0.0,0.0), (30.0, 60.0), ([0.0], [0.0])) if isa(dip, Array) dip = 90 .* rand(32) @@ -257,6 +279,28 @@ end @test err_d < 1e-7 end +@testset "JetSlantStackShiftSum3D, dip invariance" begin + theta = collect(-45:1.0:45) + phi = [0.0] + + A1 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = 0.0, azimuth = 0.0) + A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = 30 .* rand(64), azimuth = zeros(64)) + + m = zeros(domain(A1)) + for i = 1:129 + m[div(i,4)+20,i] = 1.0 + end + m[16,64] = 1 + m[32,96] = 1 + d1 = A1*m + d2 = A2*m + + err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1)) + + @show err_d + @test err_d < 1e-7 +end + @testset "JetSlantStackShiftSum3D, offsets parity" begin theta = collect(-45:5.0:45) phi = collect(0.0:45.0:135.0) @@ -308,7 +352,7 @@ end similarity = dot(d1, d2) / (norm(d1) * norm(d2)) @show similarity - @test similarity > 0.9 + @test similarity > 0.95 end @testset "JopSlantStackShiftSum vs JopSlantStackShiftSum3D, parity" begin