reimplement imROF

[johnnychen94]

Images.jl/src/algorithms.jl

Line 614 in 582e113

function imROF(img::AbstractMatrix{T}, λ::Number, iterations::Integer) where T<:NumberLike

If I'm not confusing, the current implementation uses the Algorithm 1 in

Getreuer, P. (2012). Rudin-Osher-Fatemi total variation denoising using split Bregman. Image Processing On Line, 2, 74-95.

In the same paper, there's a split-bregman implementation, which is a more robust and widely-used algorithm.

Goldstein, T., & Osher, S. (2009). The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), 323-343.

I'm going to reimplemente this imROF with split-bregman iteration @timholy . To do this, a psf2otf function is needed to accelerate the computation(PR: JuliaImages/ImageFiltering.jl#81)

[timholy]

I see you need the psf2otf (or whatever we want call it) to do this. I can probably do a "clean-room" implementation, I did not look at your code and the fundamental operation is pretty simple.

You'll need to give me a few days, though, I have urgent deadlines through Wednesday.

[johnnychen94]

Images.jl/src/algorithms.jl

Line 563 in 582e113

forwarddiffx(u::AbstractMatrix) = [u[:,2:end] u[:,end:end]] - u

@timholy There're five helper functions( forwarddiffx,forwarddiffy,backdiffx,backdiffy,div) to implement current imROF here. May I change/delete this when I do my own implementations?

As a simple test, I reimplemented my version of ROF_gradientdescent(by changing this), which is 5~10x faster than the current one.

For example, one of the changes is

function forwarddiffx(u::AbstractMatrix) # circular ∇x
    out = similar(u)
    out[:, 2:end] = u[:, 2:end] - u[:, 1:end-1]
    out[:, 1] = u[:, 1] - u[:, end]
    return out
end

Benchmark

julia> @time imdenoise(noisy_img,1e-3,500);
19.789209 seconds (36.89 k allocations: 27.350 GiB, 35.50% gc time)
julia> @time imROF(noisy_img,1e-3,500);
106.927849 seconds (1.33 G allocations: 40.441 GiB, 20.57% gc time)

BTW, do we have psnr/ssim implementation in Images? If not, I might help to do this.

[timholy]

Changing the algorithm is perfectly fine. Changing the result would require some kind of deprecation strategy. Or if you want to add circular boundary conditions as an option, we could follow the model of ImageFiltering.

[johnnychen94]

@timholy are you available these days to implement the psf2otf and otf2psf utils? 😄 I might need this in the near future.

[timholy]

Sure. Keep in mind that unless you want to change sizes, psf2otf is just called fft 😄 :

julia> using ImageFiltering

julia> kern = Kernel.gaussian(2)
OffsetArray(::Array{Float64,2}, -4:4, -4:4) with eltype Float64 with indices -4:4×-4:4:
 0.000763447  0.00183141  0.00342153  …  0.00342153  0.00183141  0.000763447
 0.00183141   0.00439334  0.00820783     0.00820783  0.00439334  0.00183141 
 0.00342153   0.00820783  0.0153342      0.0153342   0.00820783  0.00342153 
 0.0049783    0.0119423   0.0223112      0.0223112   0.0119423   0.0049783  
 0.00564116   0.0135324   0.0252819      0.0252819   0.0135324   0.00564116 
 0.0049783    0.0119423   0.0223112   …  0.0223112   0.0119423   0.0049783  
 0.00342153   0.00820783  0.0153342      0.0153342   0.00820783  0.00342153 
 0.00183141   0.00439334  0.00820783     0.00820783  0.00439334  0.00183141 
 0.000763447  0.00183141  0.00342153     0.00342153  0.00183141  0.000763447

julia> using FFTW

julia> fft(kern)
9×9 Array{Complex{Float64},2}:
         1.0+0.0im          …     0.310249-0.26033im    
    0.310249+0.26033im            0.164026-1.38778e-17im
  0.00174698+0.00990761im       0.00312125+0.00261904im 
 -0.00254613+0.00441003im      0.000358128+0.00203104im 
  0.00133217-0.000484872im     0.000287079-0.000497236im
  0.00133217+0.000484872im  …  0.000539532-0.000196374im
 -0.00254613-0.00441003im        -0.001938-0.000705374im
  0.00174698-0.00990761im      -0.00203725-0.00352862im 
    0.310249-0.26033im           0.0284829-0.161534im   

That's because we can have arbitrary axes, which is not supported in Python & Matlab.

[johnnychen94]

Thanks! psf2otf is used here to accelerate the iteration incorporating with the convolution theorem, and changing sizes is needed to do element-wise multiplication.

# split-bregman iteration, aka ADMM
∇x = psf2otf([1,-1],size(clean_img))
∇y = psf2otf([1;-1],size(clean_img))
∇T∇ = @. abs(∇x)^2 + abs(∇y)^2

d1 = 0 .* u0; b1 = d1; 
d2 = d1; b2 = d2;
u = u0;
λ = 0.05; μ = 0.05;

for iter=1:200
    Nomin = λ*u0 .+ μ .* (backdiffx(d1 .- b1) .+ backdiffy(d2 .- b2));
    Denom = λ .+ μ .* ∇T∇;
    u = real(ifft(fft(Nomin)./Denom));
    
    [d1,d2] = shrink2(Dx(u) .+ b1, Dy(u) .+ b2, 1 ./ mu);
    b1 = b1 .+ Dx(u) .- d1;
    b2 = b2 .+ Dy(u) .- d2;
end

u = real(ifft(fft(Nomin)./Denom)) is the solution to the following equation （symbols are not identical)

Goldstein, T., & Osher, S. (2009). The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), 323-343.

[timholy]

I guess a key question is whether we want to target numerical applications or physical concepts. From a conceptual standpoint, the OTF is a function of frequencies that have physical units: e.g., if your pixels have 1mm spacing then the highest frequency in your OTF is 0.5/mm. Should we return this as an AxisArray? ImageFiltering currently doesn't depend on AxisArrays, but Images does, so it could be implemented here (ref JuliaImages/ImageFiltering.jl#81).

Perhaps an alternative would be to return the AxisArray only if the input has AxisArray-axes, i.e., if the user has gone to some trouble to specify non-default pixelspacing.

[timholy]

Another question: would you rather use rfft and irfft rather than fft/ifft?

[johnnychen94]

whether we want to target numerical applications or physical concepts

Will the physical concept "Optical Transfer function" more appropriate to place into Optical packages such as QuantumOptics.jl? 🤔 In my use case, perhaps rfft(x [, dims] [, outsize])) makes more sense

ImageFiltering currently doesn't depend on AxisArrays, but Images does, so it could be implemented here

I guess ImageFiltering is more appropriate, even if we implement it here in Images, we might move it to ImageFiltering eventually since Images is quite a huge dependency. Or maybe it can be implemented in FFTW, I don't know how commonly/broadly this would be used.

Another question: would you rather use rfft and irfft rather than fft/ifft?

rfft and irfft are great! I wasn't aware of them until now. This piece of codes was written several months ago.

[timholy]

OK, let's start with a "numeric" implementation and add others as needed.