Product Operator

This operator describes the product or composition between two operators:

weights = collect(range(0, 1, length = N*N))
wop = WeightingOp(weights)
fop = FFTOp(ComplexF64, shape = (N, N));

A feature of LinearOperators.jl is that operator can be cheaply transposed, conjugated and multiplied and only in the case of a matrix-vector product the combined operation is evaluated.

tmp_op = wop * fop
tmp_freqs = tmp_op * vec(image)

65536-element Vector{ComplexF64}:
                   0.0 + 0.0im
 3.0462000080526296e-8 - 1.509759737958227e-7im
 5.2459024407156914e-8 - 1.3584434047954445e-7im
 1.9100845391789574e-8 + 3.024197548264238e-7im
 -7.303894654075166e-8 + 5.991128465798023e-7im
 -1.260072820095197e-7 + 1.3330537040721718e-7im
 -3.123155054383818e-8 - 8.755270192366595e-7im
 1.1135061314377679e-7 - 6.874887637015389e-7im
 1.2909126790827903e-7 + 1.0000116481888043e-6im
 1.3301377861841845e-7 + 1.3770225575294798e-6im
                       ⋮
  -0.03197243815052351 + 0.00948452668203438im
 -0.022453858786396375 + 0.0029033444279687043im
   0.01791066659457831 - 0.007144139531396278im
  0.031171127533136264 - 0.010672659555665206im
   0.00720835564263728 - 0.003456934510851998im
  -0.02567835205255712 + 0.006717289791803288im
  -0.02183521904611276 + 0.00494607648124139im
  0.021630426098163028 - 0.003264684881024579im
  0.040377583935354844 - 0.004040561564771873im

Similar to the WeightingOp, the main difference with the product operator provided by LinearOperatorCollection is the dedicated type, which allows for code specialisation.

pop = ProdOp(wop, fop)
typeof(pop)

ProdOp{ComplexF64, WeightingOp{Float64, Vector{Float64}}, LinearOperatorFFTWExt.FFTOpImpl{ComplexF64, Vector{ComplexF64}, FFTW.cFFTWPlan{ComplexF64, -1, true, 2, Tuple{Int64, Int64}}, FFTW.cFFTWPlan{ComplexF64, 1, true, 2, Tuple{Int64, Int64}}}, Vector{ComplexF64}}

and the ability to retrieve the components:

typeof(pop.A)

WeightingOp{Float64, Vector{Float64}}

and

typeof(pop.B)

LinearOperatorFFTWExt.FFTOpImpl{ComplexF64, Vector{ComplexF64}, FFTW.cFFTWPlan{ComplexF64, -1, true, 2, Tuple{Int64, Int64}}, FFTW.cFFTWPlan{ComplexF64, 1, true, 2, Tuple{Int64, Int64}}}

Otherwise they compute the same thing:

pop * vec(image) == tmp_op * vec(image)

true

Note that a product operator is not thread-safe.

We can again visualize our result:

weighted_frequencies = reshape(pop * vec(image), N, N)
image_inverse = reshape(adjoint(pop) * vec(weighted_frequencies), N, N)
fig = Figure()
plot_image(fig[1,1], image, title = "Image")
plot_image(fig[1,2], abs.(weighted_frequencies) .+ eps(), title = "Frequency Domain", colorscale = log10)
plot_image(fig[1,3], real.(image_inverse), title = "Image after Inverse")
resize_to_layout!(fig)
fig

This page was generated using Literate.jl.