The principles of Images.jl: Part I

This blog demonstrates how we (my mentor @zygmuntszpak and me @johnnychen94) plan to evolve the entire ecosystem in my JSoC project "Towards Better Images.jl Ecosystem". As a side product of it, a new package ImageQualityIndexes.jl is released.

In this blog, I'll make a short introduction of ImageQualityIndexes.jl, then explain how we make all the decisions.

A quick demo of ImageQualityIndexes.jl

Before explaining the coding idea in details, I'd like to show how we can use this package.

# Let's add and import some package first
pkg"activate ."
pkg"add ImageQualityIndexes ImageCore ImageDistances ImageFiltering ImageShow TestImages"
using ImageCore, ImageDistances, ImageShow, TestImages, ImageFiltering, ImageQualityIndexes, Test
# ImageQualityIndexes will be reexported by Images in the future

115.7s

Julia

Assume that you have implemented a fantastic image restoration algorithm, which in our case is mean filter, we want to show how fantastic this algorithm is. Two most famous indexes/methods used here are: Peak signal-to-noise ratio (PSNR) and Structural similarity (SSIM). They're implemented in ImageQualityIndexes.jl now🎉

# let's load lena with great respect, who significantly pushes the progress of image processing.😄
img = testimage("lena_gray_256") .|> float64
noisy_img = colorview(Gray, img .+ 0.1 .* randn(size(img))) # add gaussian noise of std 0.1

kernel = centered(ones(3, 3)./9)
denoised_img = imfilter(noisy_img, kernel) # denoise image with mean filter

# great, it works
# now let's see how fantastic it is:

# 1. quantity performance: psnr, ssim
psnrs = [psnr(img, img), psnr(noisy_img, img), psnr(denoised_img, img)]
ssims = [ssim(img, img), ssim(noisy_img, img), ssim(denoised_img, img)]
println("PSNRs:", psnrs)
println("SSIMs:", ssims)

# 2. visual quality performance
[img noisy_img denoised_img]

28.5s

Julia

As your algorithm becoming better, PSNR goes to Inf and SSIM goes to 1.

That's how we use it in most cases, quite easy, right? Actually, it's also made robust and smart enough, of course, in a simple way.

Review different image types

If you did some real image processing work in other languages, you probably know that

there is a uint8 storage type with values in {0, 1, ..., 255}, and a double storage type with values in [0, 1]. Don't forget that there're uint16, etc.
there's a double(img) which simply does a type conversion and im2double(img) that does additional value mapping from {0, 1, ..., 255} to [0,1]
there is a m*n*3 RGB image with red(R), green(G) and blue(B) channels, and a m*n*3 Lab image with lightness(L) and a, b channels. No doubt there are more.

It immediately becomes a nightmare for package developers and package users, that nobody knows what exactly a m*n*3 array represents. Is it Lab image, RGB image or even 3D gray image?

In short, as you can imagine, flexibility and robustness seems conflict with each other.

How Images.jl save your life?

Luckily, this is a lovely story in JuliaImages, thanks to the work from @timholy and all other contributors in the dark age of Julia.

In JuliaImages, we have an intermediate type Colorant that mimics the meaning of pixel; all images are of type AbstractArray{<:Colorant, N}. This notation is extremely powerful in the following senses:

we immediately know that Array{RGB, 2} is a 2D RGB image, that Array{Lab, 2} is a 2D Lab image, and that Array{Gray, 3} is a 3D Gray image.
we can easily write flexible and robust algorithms using the multiple dispatches on eltype(img).

To ease our life by avoiding unnecessary ambiguities, we made the following decisions:

Rule of thumb:

AbstractArray{<:Number} is treated as gray images by default
abandon UIntxx type in favor of FixedPointNumbers.Normed

Take the example of PSNR, the definition is quite simple for gray images

# definition of PSNR
function psnr_1(x::AbstractArray{T}, ref::AbstractArray{T}, peakval) where T <: Union{AbstractGray, Number}
  20log10(peakval) - 10log10(mse(x, ref))
end

# RGB images are treated as 3D Gray images, so we only need to unwrap it.
function psnr_1(x::AbstractArray{T}, ref::AbstractArray{T}, peakval) where T <: AbstractRGB
  psnr_1(channelview(x), channelview(ref), peakval)
end

@test psnr_1(noisy_img, img, 1) == psnr(noisy_img, img)

2.1s

Julia

Test Passed

Since we don't need to guess the input types by checking the array value, which consists of a major part of codes in other languages, we need less codes here. Also, it's as robust as you can imagine; if you pass Lab images, you will get a MethodError. No we don't guess the inputs in JuliaImages; a lot of type annotations are used in Images.jl to avoid ambiguity.

In other languages, we also need to guess the peak signal value, while in Julia, it's just an additional three lines of code:

# infer the peakvalue according to pixel type
peakval(::Type{T}) where T<:AbstractGray = 1
peakval(::Type{T}) where T<:AbstractRGB = 1
psnr_1(x, ref) = psnr_1(x, ref, peakval(eltype(ref)))

@test psnr_1(noisy_img, img) == psnr(noisy_img, img) # great!

0.6s

Julia

Test Passed

Is that all you're expecting? Good. ImageQualityIndexes is designed in this simple, flexible and robust way. We are going to make this happen to other JuliaImages packages.

Oh there's another demo

Recall how we calculate psnrs in the demo:

psnrs = [psnr(img, img), psnr(noisy_img, img), psnr(denoised_img, img)]
ssims = [ssim(img, img), ssim(noisy_img, img), ssim(denoised_img, img)]
println("PSNRs:", psnrs)
println("SSIMs:", ssims)

0.9s

Julia

One may ask: can we make it shorter? Yes we can do this using list comprehension.

psnrs = [psnr(x, img) for x in (img, noisy_img, denoised_img)]
ssims = [ssim(x, img) for x in (img, noisy_img, denoised_img)]
println("PSNRs:", psnrs)
println("SSIMs:", ssims)

1.8s

Julia

Is this the shortest? Definitely NO.

psnrs, ssims = [[assess(iqi, x, img) for x in (img, noisy_img, denoised_img)]
                                     for iqi in (PSNR(), SSIM())]
println("PSNRs:", psnrs)
println("SSIMs:", ssims)

1.6s

Julia

That said, the last one isn't as easy to be understood as the first two. But you may notice another way of calling functions in ImageQualityIndexes.jl through assess. Actually, there's one more, PSNR()(noisy_img, img), which sometimes is called a functor in other languages, i.e., a class instance that behaves like a function.

# three set of APIs in ImageQualityIndexes.jl
iqi = PSNR()
@test assess(iqi, noisy_img, img) == iqi(noisy_img, img) == psnr(noisy_img, img)

0.4s

Julia

true

Why we invent a new name?

I bet you must be wondering why we need to invent a new name `assess`? Our reason for it is:

> Readability counts.

Still take the example of ImageQualityIndexes.jl, but please do think it in a more general way:

Although assess and functor work equivalently,

assess(iqi, x, ref) reads as "assess the image quality of x using method iqi with information ref"
iqi(x, ref) is not easy to be understood at the first look, since one might guess what iqi is -- It might be a Flux model though :D

Although assess and psnr/ssim work equivalently,

From assess(SSIM(kernel, W), x, ref) one immediately knows that kernel, W are arguments to SSIM without even looking into the docsting.
From ssim(x, ref, kernel, W) we don't know this fact unless we check the docstring, if there is a well-documented one. (check ?ssim for it)

So far I think assess is the most explicit and accurate API.

Don't worry, the codes are still pretty simple even if we have three sets of APIs

abstract type ImageQualityIndex_1 end

# we use functor to link assess to implementations
assess_1(iqi::ImageQualityIndex_1, x, ref) = iqi(x, ref)

struct PSNR_1 <:ImageQualityIndex_1 end
peakval(::Type{T}) where T<:AbstractGray = 1
peakval(::Type{T}) where T<:AbstractRGB = 1

# hide implementation details in another not-exported name
(iqi::PSNR_1)(x, ref, peakval) = _psnr(x, ref, peakval)
(iqi::PSNR_1)(x, ref) = iqi(x, ref, peakval(eltype(ref)))

psnr_1(x, ref, peakval) = _psnr(x, ref, peakval)
psnr_1(x, ref) = psnr_1(x, ref, peakval(eltype(ref)))

function _psnr(x::AbstractArray{T}, ref::AbstractArray{T}, peakval) where T <: Union{AbstractGray, Number}
	20log10(peakval) - 10log10(mse(x, ref))
end
function _psnr(x::AbstractArray{T}, ref::AbstractArray{T}, peakval) where T <: AbstractRGB
  _psnr(channelview(x), channelview(ref), peakval)
end

1.0s

Julia

_psnr (generic function with 2 methods)

The following figure shows how I design this package.

@test assess_1(PSNR_1(), noisy_img, img) == assess(PSNR(), noisy_img, img)
@test PSNR_1()(noisy_img, img) == PSNR()(noisy_img, img)

1.8s

Julia

Test Passed

How we understand the three sets of APIs?

assess is the standard way, aka protocol, of calling any image quality assessment method
functors like iqi(x, ref) belongs to implementation details that links different algorithms to assess without much coupling.
convenient names psnr and ssim are made for convenient and compatibility usage. After all, users from MATLAB are more familiar to psnr and ssim.

Conclusion

This is going to be the end of the story. Not all things are organized logically, but I still hope you now understand the basic and noteworthy principles of JuliaImages. Just to make it like a conclusion, they're:

AbstractArray{<:Number} is treated as gray images by default
no UIntxx
take the advantage of dispatch, and, don't overuse dispatch
Readability counts