# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1610085490 18000
# Node ID 1b9e90ca81e3887755ef0d9666fb73796e7df4ae
# Parent  b32ced049cbd1365e27f969d4ce272e5a3080b24
Add FilterKernel LinOp

diff -r b32ced049cbd -r 1b9e90ca81e3 src/ImFilter.jl
--- a/src/ImFilter.jl	Fri Jan 08 00:34:47 2021 -0500
+++ b/src/ImFilter.jl	Fri Jan 08 00:58:10 2021 -0500
@@ -8,27 +8,32 @@
 module ImFilter
 
 using OffsetArrays
-using AlgTools.Util: @threadsif
+using AlgTools.Util: @threadsif, norm₁
+using AlgTools.LinOps
 
 ##########
 # Exports
 ##########
 
 export simple_imfilter,
-       gaussian
+       simple_imfilter!,
+       simple_imfilter_adjoint,
+       simple_imfilter_adjoint!,
+       gaussian,
+       FilterKernel
 
 ##############
 # The routine
 ##############
 
+Image = Array{Float64,2}
+Kernel = OffsetArray{Float64,2,Image}
+
 @inline function inside(i, aʹ, bʹ, a, b)
      return (max(a, i - aʹ) - i):(min(b,  i + bʹ) - i)
 end
 
-function simple_imfilter(b::Array{Float64,2},
-                         kernel::OffsetArray{Float64,2,Array{Float64,2}};
-                         threads::Bool=true)
-
+function simple_imfilter!(res::Image, b::Image, kernel::Kernel; threads::Bool=true)
     n, m = size(b)
     k, 𝓁 = size(kernel)
     o₁, o₂ = kernel.offsets
@@ -36,9 +41,7 @@
     b₁, b₂ = -1 - o₁, -1 - o₂
     kp = kernel.parent
 
-    @assert(isodd(k) && isodd(𝓁))
-
-    res = similar(b)
+    @assert(isodd(k) && isodd(𝓁) && size(res)==size(b))
 
     @threadsif threads for i=1:n
         @inbounds for j=1:m
@@ -57,6 +60,72 @@
     return res
 end
 
+function simple_imfilter(b::Image, kernel::Kernel; threads::Bool=true)
+    res = similar(b)
+    simple_imfilter!(res, b, kernel)
+end
+
+function simple_imfilter_adjoint!(res::Image, b::Image, kernel::Kernel; threads::Bool=true)
+    n, m = size(b)
+    k, 𝓁 = size(kernel)
+    o₁, o₂ = kernel.offsets
+    a₁, a₂ = k + o₁, 𝓁 + o₂
+    b₁, b₂ = -1 - o₁, -1 - o₂
+    kp = kernel.parent
+
+    @assert(isodd(k) && isodd(𝓁) && size(res)==size(b))
+
+    res .= 0
+
+    @threadsif threads for i=1:n
+        @inbounds for j=1:m
+            it₁ = inside(i, a₁, b₁, 1, n)
+            it₂ = inside(j, a₂, b₂, 1, m)
+            for p=it₁
+                @simd for q=it₂
+                    res[i+p,j+q] += kp[p-o₁, q-o₂]*b[i, j]
+                end
+            end
+        end
+    end
+
+    return res
+end
+
+function simple_imfilter_adjoint(b::Image, kernel::Kernel; threads::Bool=true)
+    res = similar(b)
+    simple_imfilter_adjoint!(res, b, kernel)
+end
+
+###########################
+# Abstract linear operator
+###########################
+
+struct FilterKernel <: AdjointableOp{Image, Image}
+    kernel::Kernel
+end
+
+function (op::FilterKernel)(b::Image)
+    return simple_imfilter(b, op.kernel)
+end
+
+function LinOps.inplace!(y::Image, op::FilterKernel, x::Image)
+    return simple_imfilter!(y, x, op.kernel)
+end
+
+function LinOps.calc_adjoint(op::FilterKernel, y::Image)
+    return simple_imfilter_adjoint(y, op.kernel)
+end
+
+function LinOps.calc_adjoint!(res::Image, op::FilterKernel, y::Image)
+    return simple_imfilter_adjoint!(res, y, op.kernel)
+end
+
+function LinOps.opnorm_estimate(op::FilterKernel)
+    # Due to |f * g|_p ≤ |f|_p|g|_1
+    return norm₁(op.kernel)
+end
+
 ######################################################
 # Distributions. Just to avoid the long load times of
 # ImageFiltering and heavy dependencies on FFTW etc.