Mercurial > repos > ImageTools / file comparison
comparison: src/ImFilter.jl
src/ImFilter.jl
- changeset 48
- 1b9e90ca81e3
- parent 32
- 41d13bf7d637
equal
deleted
inserted
replaced
| 47:b32ced049cbd | 48:1b9e90ca81e3 |
|---|---|
| 6 __precompile__() | 6 __precompile__() |
| 7 | 7 |
| 8 module ImFilter | 8 module ImFilter |
| 9 | 9 |
| 10 using OffsetArrays | 10 using OffsetArrays |
| 11 using AlgTools.Util: @threadsif | 11 using AlgTools.Util: @threadsif, norm₁ |
| 12 using AlgTools.LinOps | |
| 12 | 13 |
| 13 ########## | 14 ########## |
| 14 # Exports | 15 # Exports |
| 15 ########## | 16 ########## |
| 16 | 17 |
| 17 export simple_imfilter, | 18 export simple_imfilter, |
| 18 gaussian | 19 simple_imfilter!, |
| 20 simple_imfilter_adjoint, | |
| 21 simple_imfilter_adjoint!, | |
| 22 gaussian, | |
| 23 FilterKernel | |
| 19 | 24 |
| 20 ############## | 25 ############## |
| 21 # The routine | 26 # The routine |
| 22 ############## | 27 ############## |
| 23 | 28 |
| 29 Image = Array{Float64,2} | |
| 30 Kernel = OffsetArray{Float64,2,Image} | |
| 31 | |
| 24 @inline function inside(i, aʹ, bʹ, a, b) | 32 @inline function inside(i, aʹ, bʹ, a, b) |
| 25 return (max(a, i - aʹ) - i):(min(b, i + bʹ) - i) | 33 return (max(a, i - aʹ) - i):(min(b, i + bʹ) - i) |
| 26 end | 34 end |
| 27 | 35 |
| 28 function simple_imfilter(b::Array{Float64,2}, | 36 function simple_imfilter!(res::Image, b::Image, kernel::Kernel; threads::Bool=true) |
| 29 kernel::OffsetArray{Float64,2,Array{Float64,2}}; | |
| 30 threads::Bool=true) | |
| 31 | |
| 32 n, m = size(b) | 37 n, m = size(b) |
| 33 k, 𝓁 = size(kernel) | 38 k, 𝓁 = size(kernel) |
| 34 o₁, o₂ = kernel.offsets | 39 o₁, o₂ = kernel.offsets |
| 35 a₁, a₂ = k + o₁, 𝓁 + o₂ | 40 a₁, a₂ = k + o₁, 𝓁 + o₂ |
| 36 b₁, b₂ = -1 - o₁, -1 - o₂ | 41 b₁, b₂ = -1 - o₁, -1 - o₂ |
| 37 kp = kernel.parent | 42 kp = kernel.parent |
| 38 | 43 |
| 39 @assert(isodd(k) && isodd(𝓁)) | 44 @assert(isodd(k) && isodd(𝓁) && size(res)==size(b)) |
| 40 | |
| 41 res = similar(b) | |
| 42 | 45 |
| 43 @threadsif threads for i=1:n | 46 @threadsif threads for i=1:n |
| 44 @inbounds for j=1:m | 47 @inbounds for j=1:m |
| 45 tmp = 0.0 | 48 tmp = 0.0 |
| 46 it₁ = inside(i, a₁, b₁, 1, n) | 49 it₁ = inside(i, a₁, b₁, 1, n) |
| 53 res[i, j] = tmp | 56 res[i, j] = tmp |
| 54 end | 57 end |
| 55 end | 58 end |
| 56 | 59 |
| 57 return res | 60 return res |
| 61 end | |
| 62 | |
| 63 function simple_imfilter(b::Image, kernel::Kernel; threads::Bool=true) | |
| 64 res = similar(b) | |
| 65 simple_imfilter!(res, b, kernel) | |
| 66 end | |
| 67 | |
| 68 function simple_imfilter_adjoint!(res::Image, b::Image, kernel::Kernel; threads::Bool=true) | |
| 69 n, m = size(b) | |
| 70 k, 𝓁 = size(kernel) | |
| 71 o₁, o₂ = kernel.offsets | |
| 72 a₁, a₂ = k + o₁, 𝓁 + o₂ | |
| 73 b₁, b₂ = -1 - o₁, -1 - o₂ | |
| 74 kp = kernel.parent | |
| 75 | |
| 76 @assert(isodd(k) && isodd(𝓁) && size(res)==size(b)) | |
| 77 | |
| 78 res .= 0 | |
| 79 | |
| 80 @threadsif threads for i=1:n | |
| 81 @inbounds for j=1:m | |
| 82 it₁ = inside(i, a₁, b₁, 1, n) | |
| 83 it₂ = inside(j, a₂, b₂, 1, m) | |
| 84 for p=it₁ | |
| 85 @simd for q=it₂ | |
| 86 res[i+p,j+q] += kp[p-o₁, q-o₂]*b[i, j] | |
| 87 end | |
| 88 end | |
| 89 end | |
| 90 end | |
| 91 | |
| 92 return res | |
| 93 end | |
| 94 | |
| 95 function simple_imfilter_adjoint(b::Image, kernel::Kernel; threads::Bool=true) | |
| 96 res = similar(b) | |
| 97 simple_imfilter_adjoint!(res, b, kernel) | |
| 98 end | |
| 99 | |
| 100 ########################### | |
| 101 # Abstract linear operator | |
| 102 ########################### | |
| 103 | |
| 104 struct FilterKernel <: AdjointableOp{Image, Image} | |
| 105 kernel::Kernel | |
| 106 end | |
| 107 | |
| 108 function (op::FilterKernel)(b::Image) | |
| 109 return simple_imfilter(b, op.kernel) | |
| 110 end | |
| 111 | |
| 112 function LinOps.inplace!(y::Image, op::FilterKernel, x::Image) | |
| 113 return simple_imfilter!(y, x, op.kernel) | |
| 114 end | |
| 115 | |
| 116 function LinOps.calc_adjoint(op::FilterKernel, y::Image) | |
| 117 return simple_imfilter_adjoint(y, op.kernel) | |
| 118 end | |
| 119 | |
| 120 function LinOps.calc_adjoint!(res::Image, op::FilterKernel, y::Image) | |
| 121 return simple_imfilter_adjoint!(res, y, op.kernel) | |
| 122 end | |
| 123 | |
| 124 function LinOps.opnorm_estimate(op::FilterKernel) | |
| 125 # Due to |f * g|_p ≤ |f|_p|g|_1 | |
| 126 return norm₁(op.kernel) | |
| 58 end | 127 end |
| 59 | 128 |
| 60 ###################################################### | 129 ###################################################### |
| 61 # Distributions. Just to avoid the long load times of | 130 # Distributions. Just to avoid the long load times of |
| 62 # ImageFiltering and heavy dependencies on FFTW etc. | 131 # ImageFiltering and heavy dependencies on FFTW etc. |
