Image samples and derivatives for direct/photometric image alignment

Specific Domains Signal and Image Processing
1

I’m planning to do some direct (or photometric) image alignment (or registration) in Julia. By this I mean methods based on the Lucas-Kanade method. For this I’ll need not only image samples at non-integer locations, but also the derivatives of these values. Has anyone implemented any optimizations that use both image values and derivatives? If so, how did you do it? Are there some image samplers/interpolators that provide derivatives, or overload calls with autodiff types (e.g. ForwardDiff dual variables), or some Chain Rules?

2

Interpolations.jl provides gradients out of the box, including for RGB pixel types:

julia> using ImageCore, Interpolations

julia> img = [RGB(0,0,0) RGB(0,1,0); RGB(0,0,1) RGB(1,1,1)]
2×2 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(0.0,0.0,0.0)  RGB{N0f8}(0.0,1.0,0.0)
 RGB{N0f8}(0.0,0.0,1.0)  RGB{N0f8}(1.0,1.0,1.0)

julia> itp = linear_interpolation(axes(img), img);

julia> Interpolations.gradient(itp, 1.5, 1.5)
2-element StaticArraysCore.SVector{2, RGB{Float64}} with indices SOneTo(2):
 RGB{Float64}(0.5,0.0,1.0)
 RGB{Float64}(0.5,1.0,0.0)
2 Likes
3

Thank you! Strangely this function doesn’t appear in the Public API page of Interpolations’ docs. Also, it seems like you can’t get values and gradients in the same call, which is a shame, since much of the computation is duplicated.

4

Gradients for images are available in ImageFiltering.jl too with imgradients and this might be of interest too: https://github.com/JuliaImages/ImageBase.jl/blob/master/src/diff.jl .

5

Thank you. However, for getting samples and gradients at non-integer locations, I don’t believe computing gradient images first is a good idea, for the following reasons:

  • There is a non-negligible computational overhead to first computing gradient images. If you only need a few samples, then this will be swamped by the pre-computation.
  • Computing gradient images (then sampling these at non-integer locations) actually triples the amount of data you need to hold in memory, and also almost triples the amount of work that needs to be done when computing values and gradients at non-integer locations. Since you’re also pulling 3 times more data through the cache, it also leads to more cache misses, slowing things down further.
  • Gradients computed this way will not be analytic derivatives of the sampled values, but rather some approximation. This means that an optimization that uses these values is not guaranteed to converge, or may converge more slowly.
1 Like
6

BTW I was seeing the docs of ImageCore.jl today and these traits might be useful to you: Traits · ImageCore

1 Like
7

The computation itself is cheap enough that I’d expect it to be dominated by memory bandwidth if you calculate both the value and gradient at the same time

julia> function gradient_only(itp_img, locs)
           s = itp_img(locs[1]...)*0
           @inbounds for loc in locs
               s += sum(Interpolations.gradient(itp_img, loc...))
           end
           return s
       end
gradient_only (generic function with 1 method)

julia> function gradient_value(itp_img, locs)
           s = itp_img(locs[1]...)*0
           @inbounds for loc in locs
               s += itp_img(loc...)
               s += sum(Interpolations.gradient(itp_img, loc...))
           end
           return s
       end
gradient_value (generic function with 1 method)

julia> img = rand(RGB{N0f8}, 1024, 1024);

julia> locs = collect(zip(rand(axes(img, 1), 2048), rand(axes(img, 2), 2048)));

julia> itp_img = linear_interpolation(axes(img), img);

julia> @btime gradient_only($itp_img, $locs)
  69.400 μs (0 allocations: 0 bytes)
RGB{Float64}(-1.1529411764705895,-18.10980392156866,-1.8392156862744766)

julia> @btime gradient_value($itp_img, $locs)
  65.600 μs (0 allocations: 0 bytes)
RGB{Float64}(1030.321568627451,1024.2039215686275,1020.3254901960775)

Strange - calculating the value and gradient was slightly faster than calculating the gradient alone! Probably a benchmarking artifact.