GitHub

ImageGeoms overview

This page explains the Julia package ImageGeoms.

This page comes from a single Julia file: 1-overview.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: 1-overview.ipynb, or open it in binder here: 1-overview.ipynb.

Setup

Packages needed here.

using ImageGeoms: ImageGeom, axis, axisf, axesf, downsample, grids, oversample
import ImageGeoms # ellipse
using AxisArrays
using MIRTjim: jim, prompt
using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)
using Unitful: mm, s
using InteractiveUtils: versioninfo

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

isinteractive() ? jim(:prompt, true) : prompt(:draw);

Overview

When performing tomographic image reconstruction, one must specify the geometry of the grid of image pixels. (In contrast, for image denoising and image deblurring problems, one works with the given discrete image and no physical coordinates are needed.)

The key parameters of a grid of image pixels are

  • the size (dimensions) of the grid, e.g., 128 × 128,
  • the spacing of the pixels, e.g., 1mm × 1mm,
  • the offset of the pixels relative to the origin, e.g., (0,0)

The data type ImageGeom describes such a geometry, for arbitrary dimensions (2D, 3D, etc.).

There are several ways to construct this structure. The default is a 128 × 128 grid with pixel size $\Delta_X = \Delta_Y = 1$ (unitless) and zero offset:

ig = ImageGeom()
ImageGeom{2, NTuple{2,Float32}, FillArrays.Trues{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{2,Int64} (128, 128)
 deltas::NTuple{2,Float32} (1.0f0, 1.0f0)
 offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
 mask: 128×128 Ones{Bool} {16384 of 16384}

Here is a 3D example with non-cubic voxel size:

ig = ImageGeom( (512,512,128), (1,1,2), (0,0,0) )
ImageGeom{3, NTuple{3,Int64}, FillArrays.Trues{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{3,Int64} (512, 512, 128)
 deltas::NTuple{3,Int64} (1, 1, 2)
 offsets::NTuple{3,Float32} (0.0f0, 0.0f0, 0.0f0)
 mask: 512×512×128 Ones{Bool} {33554432 of 33554432}

To avoid remembering the order of the arguments, named keyword pairs are also supported:

ig = ImageGeom( dims=(512,512,128), deltas=(1,1,2), offsets=(0,0,0) )
ImageGeom{3, NTuple{3,Int64}, FillArrays.Trues{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{3,Int64} (512, 512, 128)
 deltas::NTuple{3,Int64} (1, 1, 2)
 offsets::NTuple{3,Float32} (0.0f0, 0.0f0, 0.0f0)
 mask: 512×512×128 Ones{Bool} {33554432 of 33554432}

Units

The pixel dimensions deltas can (and should!) be values with units.

Here is an example for a video (2D+time) with 12 frames per second:

ig = ImageGeom( dims=(640,480,1000), deltas=(1mm,1mm,(1//12)s) )
ImageGeom{3, Tuple{Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Unitful.Quantity{Rational{Int64}, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, FillArrays.Trues{3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{3,Int64} (640, 480, 1000)
 deltas::Tuple{Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Unitful.Quantity{Rational{Int64}, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}} (1 mm, 1 mm, 1//12 s)
 offsets::NTuple{3,Float32} (0.0f0, 0.0f0, 0.0f0)
 mask: 640×480×1000 Ones{Bool} {307200000 of 307200000}

Methods

An ImageGeom object has quite a few methods; axes and axis are especially useful:

ig = ImageGeom( dims=(7,8), deltas=(3,2), offsets=(0,0.5) )
axis(ig, 2)
-8.0f0:2.0f0:6.0f0

or an axis of length n with spacing Δ (possibly with units) and (always unitless but possibly non-integer) offset the axis is a subtype of AbstractRange of the form ( (0:n-1) .- ((n - 1)/2 + offset) ) * Δ

These axes are useful for plotting:

ig = ImageGeom( dims=(10,8), deltas=(1mm,1mm), offsets=(0.5,0.5) )
ImageGeom{2, NTuple{2,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, FillArrays.Trues{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{2,Int64} (10, 8)
 deltas::NTuple{2,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (1 mm, 1 mm)
 offsets::NTuple{2,Float32} (0.5f0, 0.5f0)
 mask: 10×8 Ones{Bool} {80 of 80}

_ticks(x, off) = [x[1]; # hopefully helpful tick marks
    iseven(length(x)) && iszero(off) ?
        oneunit(eltype(x)) * [-0.5,0.5] : zero(eltype(x)); x[end]]
showgrid = (ig) -> begin # x,y grid locations of pixel centers
    x = axis(ig, 1)
    y = axis(ig, 2)
    (xg, yg) = grids(ig)
    scatter(xg, yg; label="", xlabel="x", ylabel="y",
        xlims = extrema(x) .+ (ig.deltas[1] * 0.5) .* (-1,1),
        xticks = _ticks(x, ig.offsets[1]),
        ylims = extrema(y) .+ (ig.deltas[2] * 0.5) .* (-1,1),
        yticks = _ticks(y, ig.offsets[2]),
        widen = true, aspect_ratio = 1, title = "offsets $(ig.offsets)")
end
showgrid(ig)
Example block output

Axes labels can have units

prompt()

Offsets (unitless translation of grid)

The default offsets are zeros, corresponding to symmetric sampling around origin:

ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm) )
p = showgrid(ig)
Example block output
prompt();

That default for offsets is natural for tomography when considering finite pixel size:

square = (x,y,Δ,p) -> plot!(p, label="", color=:black,
    x .+ Δ[1] * ([0,1,1,0,0] .- 0.5),
    y .+ Δ[2] * ([0,0,1,1,0] .- 0.5),
)
square2 = (x,y) -> square(x, y, ig.deltas, p)
square2.(grids(ig)...)
plot!(p)
Example block output
prompt();

In that default geometry, the center (0,0) of the image is at a corner of the middle 4 pixels (for even image sizes). That default is typical for tomographic imaging (e.g., CT, PET, SPECT). One must be careful when using operations like imrotate or fft.

Odd dimensions

If an image axis has an odd dimension, then each middle pixel along that axis is centered at 0, for the default offset=0.

igo = ImageGeom( dims=(7,6) )
po = showgrid(igo)
square2 = (x,y) -> square(x, y, igo.deltas, po)
square2.(grids(igo)...)
plot!(po)
Example block output
prompt();

AxisArrays

There is a natural connection between ImageGeom and AxisArrays. Note the automatic labeling of units (when relevant) on all axes by MIRTjim.jim.

ig = ImageGeom( dims=(60,48), deltas=(1.5mm,1mm) )
za = AxisArray( ImageGeoms.ellipse(ig) * 10/mm ; x=axis(ig,1), y=axis(ig,2) )
jim(za, "AxisArray example")
Example block output
prompt();

Resizing

Often we have a target grid in mind but want coarser sampling for debugging. The downsample method is useful for this.

ig = ImageGeom( dims = (512,512), deltas = (500mm,500mm) ./ 512 )
ig_down = downsample(ig, 4)
ImageGeom{2, NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, FillArrays.Trues{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{2,Int64} (128, 128)
 deltas::NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (3.90625 mm, 3.90625 mm)
 offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
 mask: 128×128 Ones{Bool} {16384 of 16384}

Other times we want to avoid an "inverse crime" by using finer sampling to simulate data; use oversample for this.

ig_over = oversample(ig, (2,2))
ImageGeom{2, NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, FillArrays.Trues{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}
 dims::NTuple{2,Int64} (1024, 1024)
 deltas::NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (0.48828125 mm, 0.48828125 mm)
 offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
 mask: 1024×1024 Ones{Bool} {1048576 of 1048576}

Frequency domain axes

For related packages like ImagePhantoms, there are also axisf and axesf methods that return the frequency axes associated with an FFT-based approximation to a Fourier transform, where $Δ_f Δ_X = 1/N.$

ig = ImageGeom( dims=(4,5), deltas=(1mm,1mm) )
axesf(ig)
((-0.5:0.25:0.25) mm^-1, (-0.4:0.2:0.4) mm^-1)
axisf(ig, 1)
(-0.5:0.25:0.25) mm^-1
axisf(ig, 2)
(-0.4:0.2:0.4) mm^-1

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.10.3"
 "Commit 0b4590a5507 (2024-04-30 10:59 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"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/ImageGeoms.jl/ImageGeoms.jl/docs/Project.toml`
  [39de3d68] AxisArrays v0.4.7
  [e30172f5] Documenter v1.4.1
  [9ee76f2b] ImageGeoms v0.11.1 `~/work/ImageGeoms.jl/ImageGeoms.jl`
  [98b081ad] Literate v2.18.0
  [170b2178] MIRTjim v0.25.0
  [6fe1bfb0] OffsetArrays v1.14.0
  [91a5bcdd] Plots v1.40.4
  [1986cc42] Unitful v1.20.0
  [b77e0a4c] InteractiveUtils

This page was generated using Literate.jl.