GitHub

Image alignment by rank-1 method

This example illustrates 2D and 3D image alignment (estimating a simple translation between an image pair) using a rank-1 approximation to the normalized cross power spectrum, following the method of Hoge, IEEE T-MI, 2003, using the Julia language.

This page comes from a single Julia file: align1.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: align1.ipynb, or open it in binder here: align1.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "ImagePhantoms"
        "InteractiveUtils"
        "LaTeXStrings"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
        "Statistics"
        "Unitful"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using InteractiveUtils: versioninfo
using ImageGeoms: ImageGeom, axesf
using ImagePhantoms: SheppLoganEmis, spectrum, phantom
using ImagePhantoms: ellipse, ellipse_parameters
using ImagePhantoms: ellipsoid, ellipsoid_parameters
using LaTeXStrings
using LinearAlgebra: norm, svd
using MIRTjim: jim, prompt
using Plots: default, gui, plot, plot!, scatter, scatter!, savefig
using Random: seed!
using Unitful: cm, @u_str # use of physical units (cm here)
default(); default(label="", markerstrokecolor=:auto,
    guidefontsize=14, legendfontsize=14, tickfontsize=12)

The following line is helpful when running this jl-file as a script; this way it will prompt user to hit a key after each image is displayed.

isinteractive() && prompt(:prompt);

Generate data

FOV = 256cm # physical units
Nx, Ny = 128, 126
Δx = FOV / Nx # pixel size
Δy = FOV / Ny
x = ((-Nx÷2):(Nx÷2-1)) * Δx
y = ((-Ny÷2):(Ny÷2-1)) * Δy
νx = ((-Nx÷2):(Nx÷2-1)) / Nx / Δx
νy = ((-Ny÷2):(Ny÷2-1)) / Ny / Δy;

Ellipse parameters for 1st image:

param1 = ellipse_parameters(SheppLoganEmis(), fovs=(FOV,FOV))
shift2 = (1.7, 9.2) .* oneunit(Δx) # true non-integer shift
(1.7 cm, 9.2 cm)

Ellipse parameters for 2nd image:

param2 = [((p[1:2] .+ shift2)..., p[3:end]...) for p in param1];

Phantom images:

obj1 = ellipse(param1)
obj2 = ellipse(param2)
image1 = phantom(x, y, obj1)
image2 = phantom(x, y, obj2)
pim = jim(
 jim(x, y, image1, "image1"),
 jim(x, y, image2, "image2"),
 jim(x, y, image2-image1, "image2-image1")
)
Example block output
prompt()

Analytical spectra of these images (cf MRI)

spec1 = spectrum(νx, νy, obj1)
spec2 = spectrum(νx, νy, obj2);

Add noise

seed!(0)
snr2sigma(db, y) = 10^(-db/20) * norm(y) / sqrt(length(y))
σnoise = snr2sigma(40, spec1)
addnoise(y, σ) = y + σ * randn(ComplexF32, size(y))
data1 = addnoise(spec1, σnoise)
data2 = addnoise(spec2, σnoise);

Normalized cross power spectrum

ncps = @. data1 * conj(data2) / (abs(data1 * data2));

Show spectra and noisy phase difference

fun = data -> log10.(abs.(data) / maximum(abs, data1))
psp = jim(
 jim(νx, νy, fun(data1), "|spectrum1|"),
 jim(νx, νy, fun(data2), "|spectrum2|"),
 jim(νx, νy, angle.(ncps); color = :hsv, title="phase difference")
)
Example block output
prompt()

SVD

By the shift property of the 2D Fourier transform, the normalized cross power spectrum (NCPS) is

\[e^{ı 2π (ν_x d_x + ν_y d_y)} = e^{ı 2π ν_x d_x} e^{ı 2π ν_y d_y}\]

where $(d_x, d_y)$ is the 2D translation.

So the 2D NCPS (in the absence of noise) is an outer product of 1D vectors, each of which has the phase associated with the translation in one direction.

U, s, V = svd(ncps)
psig = scatter(s, xlabel=L"k", ylabel=L"σ_k", title="Scree plot")
Example block output
prompt()

Phase of principal components

u = U[:,1]
v = conj(V[:,1]) # need conjugate here because of σ u v'
puv = plot(
 plot(νx, angle.(u), xlabel=L"ν_x"),
 plot(νy, angle.(v), xlabel=L"ν_y"),
)
Example block output

We could unwrap the phase and then fit a line as suggested in the original paper.

Instead, we use a phase slope estimate described in Feiweier 2013 US Patent 8497681B2 that avoids any need for phase unwrapping.

function phase_slope(x::AbstractVector, Δν::Number; weights = 1)
   tmp = x[2:end] .* conj(x[1:end-1])
   return angle(sum(weights .* tmp)) / (2π * Δν)
end;

myshift2 = phase_slope.((u,v), 1/FOV)
_round(x) = round(u"cm", x; digits=4)
_round.(myshift2)
(1.6991 cm, 9.1996 cm)

Error: the estimated shift is remarkably close to the true shift.

error2 = myshift2 .- shift2
_round.(error2)
(-0.0009 cm, -0.0004 cm)

3D case

Here is an illustration of an extension of the method to 3D image registration.

Ellipsoid parameters for 1st image volume:

fovs = (24cm, 24cm, 20cm)
param1 = ellipsoid_parameters( ; fovs)
ob1 = ellipsoid(param1); # Vector of Ellipsoid objects

Ellipsoid parameters for 2nd image volume:

shift3 = (1.1, 2.2, 3.3) .* oneunit.(fovs)
param2 = [((p[1:3] .+ shift3)..., p[4:end]...) for p in param1];
ob2 = ellipsoid(param2);

Visualize

dims = (128, 130, 30)
ig = ImageGeom( ; dims, deltas = fovs ./ dims )
oversample = 3
image1 = phantom(axes(ig)..., ob1, oversample)
image2 = phantom(axes(ig)..., ob2, oversample)
clim = (0.95, 1.05)
plot(
 jim(axes(ig)[1:2]..., image2; title = "3D Shepp-Logan phantom slices", clim),
 jim(axes(ig)[1:2]..., image2-image1; title = "Difference", clim),
)
Example block output
prompt()

Spectra

spectrum1 = spectrum(axesf(ig)..., ob1)
spectrum2 = spectrum(axesf(ig)..., ob2);

σnoise = snr2sigma(40, spectrum1)
data1 = addnoise(spectrum1, σnoise)
data2 = addnoise(spectrum2, σnoise);

Normalized cross power spectrum

ncps = @. data1 * conj(data2) / (abs(data1 * data2));

Show spectra and noisy phase difference

fun = data -> log10.(abs.(data / maximum(abs, data1)))
iz = dims[3]÷2 .+ (0:1)
psp = jim(
 jim(axesf(ig)[1:2]..., fun(data1)[:,:,iz], "|spectrum1|"),
 jim(axesf(ig)[1:2]..., angle.(ncps[:,:,iz]), color = :hsv, title="phase difference"),
)
Example block output
prompt()

SVD for 3 different foldings of the 3D NCPS

This follows the folded NCPS method of Hoge et al., ICIP 2003, described therein in terms of the dominant singular vectors of a high-order SVD of the tensor data.

fold1 = reshape(ncps, dims[1], :)
U, s, V = svd(fold1)
u1 = U[:,1]
psig1 = scatter(s, xlabel=L"k", ylabel=L"σ_k", title="Scree plot")
pa1 = plot(angle.(u1));

fold2 = reshape(permutedims(ncps, [2 1 3]), dims[2], :)
U, s, V = svd(fold2)
u2 = U[:,1]
psig2 = scatter(s, xlabel=L"k", ylabel=L"σ_k", title="Scree plot")
pa2 = plot(angle.(u2));

fold3 = reshape(permutedims(ncps, [3 1 2]), dims[3], :)
U, s, V = svd(fold3)
u3 = U[:,1]
psig3 = scatter(s, xlabel=L"k", ylabel=L"σ_k", title="Scree plot")
pa3 = plot(angle.(u3));

p3 = plot(
 psig1, psig2, psig3,
 pa1, pa2, pa3,
 layout = (2, 3),
)
Example block output
prompt()

Estimate translation

Δν = map(i -> diff(axesf(ig)[i])[1], 1:3)
myshift3 = phase_slope.((u1,u2,u3), Δν)
_round.(myshift3)
3-element Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}:
 1.0999 cm
 2.2004 cm
 3.3001 cm

Error is small:

error3 = myshift3 .- shift3
_round.(error3)
3-element Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}:
 -0.0001 cm
  0.0004 cm
  0.0001 cm

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
 "Julia Version 1.12.5"
 "Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)"
 "Build Info:"
 "  Official https://julialang.org release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)"
 "  GC: Built with stock GC"
 "Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.6.3
  [aaaa29a8] Clustering v0.15.8
  [35d6a980] ColorSchemes v3.31.0
  [3da002f7] ColorTypes v0.12.1
  [c3611d14] ColorVectorSpace v0.11.0
  [717857b8] DSP v0.8.4
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.16.1
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.10.0
  [587475ba] Flux v0.16.9
  [a09fc81d] ImageCore v0.10.5
  [9ee76f2b] ImageGeoms v0.11.2
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.1.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.21.0
  [7035ae7a] MIRT v0.18.3
  [170b2178] MIRTjim v0.26.0
  [eb30cadb] MLDatasets v0.7.20
  [efe261a4] NFFT v0.14.3
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.10
  [429524aa] Optim v2.0.1
  [91a5bcdd] Plots v1.41.5
  [f27b6e38] Polynomials v4.1.0
  [2913bbd2] StatsBase v0.34.10
  [1986cc42] Unitful v1.28.0
  [d6d074c3] VideoIO v1.4.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0

This page was generated using Literate.jl.