JuliaImages / ImageBase.jl

MIT License
4 stars 2 forks source link

Can `restrict` be generalized to n-fold? #30

Open johnnychen94 opened 3 years ago

johnnychen94 commented 3 years ago

I'm not sure if this is possible, but I'm imagining a more efficient version of imresize(img; ratio=1//n)

timholy commented 3 years ago

In principle, yes. You'll like this, I think: restriction is the adjoint of prolongation, and prolongation is just interpolation. (https://zingale.github.io/comp_astro_tutorial/elliptic_multigrid/multigrid/multigrid-restrict-prolong.html seems vaguely useful.) So write out the prolongation as an operator and then take the adjoint.

Example:

ϕₕ(1) = ϕ[1]         # note: "indexing" with (x) is meant to convey physical coordinates different from array indexes
ϕₕ(1.5) = (ϕ[1] + ϕ[2])/2
ϕₕ(2) = ϕ[2]

so

ϕₕ = P * ϕ

where

P = [1 0; 0.5 0.5; 0 1]

Now for restriction, just take the adjoint of P:

ϕᵣ = P'*ϕₕ

so that

ϕᵣ[1] = ϕₕ[1] + ϕₕ[2]/2
ϕᵣ[2] = ϕₕ[2]/2 + ϕₕ[3]

Note that

julia> P'*P
2×2 Matrix{Float64}:
 1.25  0.25
 0.25  1.25

which, up to a scalar multiple (divide by 1.5), is "close to" the identity matrix. (Kind of a blurred version of it.)

For n-fold prolongation, we could just linearly interpolate at ϕₕ(1+1/n), ϕₕ(1+2/n), etc. Its adjoint is (up to a scalar multiple) how one should define n-fold restriction.