Chromatic Aberration Recovery
on Arbitrary Images
by
Daniel J Blueman

A dissertation submitted to the University of Bristol
in accordance with the requirements of the degree
of Master of Research in the Faculty of Engineering

September 2011

Department of Computer Science

### Abstract

Digital imaging sensor technology has continued to outpace development in
optical technology in modern imaging systems. The resulting quality loss
attributable to lateral chromatic aberration is becoming increasingly
significant as sensor resolution increases; other classes of aberration are
less significant with classical image enhancement (e.g. sharpening), whereas
lateral chromatic aberration becomes more significant. The goals of higher-
performance and lighter lens systems drive a recent need to find new ways to
overcome resulting image quality limitations.

This work demonstrates the robust and automatic minimisation of lateral
chromatic aberration, recovering the loss of image quality using both
artificial and real-world images. A series of test images are used to validate
the functioning of the algorithm, and changes across a series of real-world
images are used to evaluate the performance of the approach.

The primary contribution of this work is introduced: a novel algorithm to
robustly minimise lateral chromatic aberration in both calibration and real-
world images. This is broken down into discrete steps and detailed. The basis
of the algorithm uses chromatic correspondences to converge a set of
distortion coefficients. Other contributions are then introduced: a second
algorithm is developed to allow correction information to be correlated to the
lens model and parameters. This information is subsequently stored in a
database to allow offline correction of unseen images. Finally, an algorithm
is developed to measure image fidelity in a way more relevant to how the Human
Visual System processes information via spatial frequency analysis.

Algorithm validation is conducted through a series of steps to ensure
correctness, then artificial and test images are analysed, and the
quantification algorithm is applied to measure improvement. Lastly, the
performance of this system is compared against prevailing methods and
analysed.

Keywords: chromatic aberration, aberrations, lens systems, image fidelity

### Acknowledgements

This work was enabled by the guidance of my supervisor, Professor Majid
Mirmehdi, whom I extend my deepest gratitude to.

Additionally, I would like to thank all my friends and colleagues at the
University of Bristol for their guidance and support in the Computer Graphics
and Computer Vision labs, including Timo Kunkel. Thanks also to Sritrakool
Waeladee for useful feedback from another discipline.

Finally, I am forever grateful to my family and friends for their moral
support.

### Declaration

I declare that the work in this dissertation was carried out in accordance
with the requirements of the University’s Regulations and Code of Practice for
Research Degree Programmes and that it has not been submitted for any other
academic award. Except where indicated by specific reference in the text, the
work is the candidate’s own work. Work done in collaboration with or with the
assistance of others is indicated as such. Any views expressed in the
dissertation are those of the author.
Daniel J Blueman

### License

This thesis and associated software is licensed under the Creative Commons
Attribution-NonCommercial-ShareAlike 3.0 Unported (CC BY-NC-SA 3.0) public
license, available at creativecommons.org.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Chromatic Aberration
    3. 1.3 Lateral Chromatic Aberration
    4. 1.4 Central argument and research objective
    5. 1.5 Contributions
    6. 1.6 Overview of methodology
    7. 1.7 Thesis layout
  2. 2 Background
    1. 2.1 Introduction
    2. 2.2 Optics
      1. 2.2.1 Paraxial optics
      2. 2.2.2 Lens systems
      3. 2.2.3 Lens parameters
      4. 2.2.4 Optical properties
      5. 2.2.5 Sensor properties
      6. 2.2.6 Sensor blooming
    3. 2.3 Image EXIF data and parameters
    4. 2.4 Lens variance
    5. 2.5 Correction overview
    6. 2.6 Test images
    7. 2.7 Literature review
      1. 2.7.1 CA and non-lens systems
      2. 2.7.2 Lens-system optical correction
      3. 2.7.3 Computational correction
    8. 2.8 Existing approaches
    9. 2.9 Summary
  3. 3 Proposed method
    1. 3.1 Introduction
    2. 3.2 Approach
      1. 3.2.1 Error function
      2. 3.2.2 Decentering calibration
    3. 3.3 Coefficient recovery
    4. 3.4 Image correction
    5. 3.5 Novel approach to correction quantification
    6. 3.6 Minimiser selection
    7. 3.7 Summary
  4. 4 Methodological evaluation
    1. 4.1 Validation
      1. 4.1.1 Chequerboard test image
      2. 4.1.2 Ray-traced scene
    2. 4.2 Performance analysis
      1. 4.2.1 Real-world test image
      2. 4.2.2 Run-time cost
      3. 4.2.3 Suite of images
    3. 4.3 Comparison with existing methods
    4. 4.4 Results discussion
    5. 4.5 Summary
  5. 5 Conclusion and further work
    1. 5.1 Review
      1. 5.1.1 Mechanics
      2. 5.1.2 Results
      3. 5.1.3 Summary
      4. 5.1.4 Limitations
    2. 5.2 Claims and contributions
    3. 5.3 Future work
      1. 5.3.1 Enhanced interpolation
      2. 5.3.2 Parameter modelling
      3. 5.3.3 Implementation performance
    4. 5.4 Summary and concluding remarks

###### List of Figures

  1. 1.1 Axial Chromatic Aberration
  2. 1.2 Axial and Lateral Chromatic Aberration
  3. 1.3 Lateral Chromatic Aberration
  4. 1.4 Cropped sections exhibiting LCA
  5. 1.5 Image correction method overview
  6. 2.1 Paraxial approximation
  7. 2.2 Dual thin lens (“doublet”) with CA correction
  8. 2.3 Snell’s law
  9. 2.4 Bayer demosaicing principle
  10. 2.5 CCD blooming principle
  11. 2.6 Real-world test images
  12. 3.1 Plane subtraction first observed by Luhmann
  13. 3.2 Cropped sections from uncorrected difference map
  14. 3.3 Coefficient recovery steps
  15. 3.4 Image correction with model database
  16. 3.5 Error varying the first and second correction coefficients
  17. 4.1 Cropped sections from test image
  18. 4.2 Change in spatial frequency in pattern images
  19. 4.3 Rendered test image
  20. 4.4 Cropped sections from rendering
  21. 4.5 Change in spatial frequency in rendered images
  22. 4.6 Cropped sections from photograph
  23. 4.7 Change in spatial frequency due to LCA correction
  24. 4.8 Runtime breakdown of algorithms
  25. 4.9 Difference map area reduction
  26. 4.10 Cropped section comparison from difference maps
  27. 4.11 Change in spatial frequency for test images
  28. 4.12 Cropped sections comparing correction methods
  29. 4.13 Change in spatial frequency in the first selected image
  30. 4.14 Cropped sections comparing correction methods
  31. 4.15 Change in spatial frequency in the second selected image
  32. 4.16 Cropped sections comparing correction methods
  33. 4.17 Change in spatial frequency in the third selected image
  34. 5.1 SIFT used to detect feature points per channel

### List of Algorithms

loa

### Abbreviations

ACA

Axial Chromatic Aberration

AHD

Adaptive Homogeneity-Directed demosaicing algorithm

BFGS

Broyden-Fletcher-Goldfarb-Shanno

CA

Chromatic Aberration

CCD

Charge Coupled Device imaging sensor

CFA

Colour Filter Array

CIE

Commission Internationale de l’Eclairage (International Lighting Commission)

CMOS

Complimentary Metal Oxide Semiconductor

DFT

Discrete Fourier Transform

DMD

Digital Micromirror Device

EXIF

EXchangeable Image File-format

FPGA

Field Programmable Gate Array

GPU

Graphics Processing Unit

HVS

Human Visual System

L-BFGS-B

Limited-memory Broyden-Fletcher-Goldfarb-Shanno-Bounded

LCA

Lateral Chromatic Aberration

PPG

Patterned Pixel Grouping demosaicing algorithm

RGB

Red-Green-Blue

SIFT

Scale-Invariant Feature Transform

SLM

Spatial Light Modulator

SLR

Single Lens Reflex camera system

UD

Ultra-low Dispersion

VNG

Variable Number of Gradients demosaicing algorithm

##  Chapter 1 Introduction

###  1.1 Motivation

Consumers, hobbyists and professionals are pushing photography into areas that
were traditionally inaccessible. The expectation and latent demand for wide-
aperture and wide-angle lenses, while retaining high-fidelity is at previously
unprecedented levels.

Modern lens design devolves to a series of cost, size and performance
tradeoffs, nominally leading to an intermediate balance of these factors.
Fixed focal length prime lenses have movable elements only for focusing,
however optical and design complexity increases substantially with variable
focal length zoom lenses which are commonplace in consumer compact cameras.
Before 2005, none of the high-zoom-ratio lenses available today were
developed. Primarily for convenience purposes, these are used to replace a
selection of prime/small-ratio lenses; a significant reduction to fidelity is
commonly lost for this convenience.

There is a continual drive to reduce costs and increase margin; reduced
complexity, material costs and weight are clear economic competitive
advantages. Due to this, there has been a growing trend to move hardware
camera functionality to in-camera software instead, reducing cost, increasing
flexibility and marketable features. The most significant example of this
would be contrast-based auto-focus; this is computed as a function over the
image frame and allows mechanical focus seeking to maximise the subject’s
captured spatial frequency. Without this, specialised auto-focus sensor chips
are needed along the optical path, increasing system complexity and costs, and
attenuating available light arriving at the image sensor, e.g. due to mirrors
or prisms.

Software implementations clearly offer far more flexibility. For example,
zone-based or face auto-focus and exposure are readily possible; these were
traditionally infeasible due to the minimal processing power inside the
camera. With these restrictions lifting, the wide scope for creating a more
flexible and optimised imaging platform able to fully exploit both hardware
and software is available.

As sensor resolution has been rapidly increasing, the optical performance has
only marginally increased, thus detail attenuated in the lens system and thus
the efficiency of the whole image-capture system is reduced, often
significantly at more extreme parameters.

Lens performance is constrained by a number of well-understood areas of
optical distortion [50]. Presently, barrel (geometrical) distortion is well
understood and has robust solutions available; chromatic aberration previously
manifested at the sub-pixel level — consequently, there is far less literature
and evolution of solutions. As sensor resolution increases, the relative
impact of chromatic aberration also increases. In general, other types of
aberration pose lower limitations to image quality and more importantly are
less sensitive to sensor resolution. Consequently, this thesis will consider
the problem of quality loss resulting from chromatic aberrations and will set
out to address this area of loss.

In this work, the term aberration is generally defined as the departure from a
theoretically ideal lens system. By extension of this, image fidelity is
defined as an image with the minimum loss due to aberration, that is free of
aberration to the fullest extent. Monochromatic is defined as features
existing at all visible wavelengths equally, such that they are invariant of
chromaticity.

From the specific areas covered, there is a clear deficit that of no current
solution which provides automatic, robust and adaptive correction for the
penalty lateral chromatic aberration imposes on general lenses; this research
aims to address this gap.

###  1.2 Chromatic Aberration

Loss in image fidelity is due to a known set of aberrations, classified as
either a monochromatic or chromatic type of aberration. They occur due to the
physical interaction of light with materials, lens design constraints and
manufacture limitations. Von Seidel rigorously defined the monochromatic types
of aberration to be astigmatic, comatic, field curvature, distortion and
spherical [95]. The chromatic type of aberration exists independently from
monochromatic aberration, and is broken down into two further types: Axial or
longitudinal Chromatic Aberration (ACA) and Lateral or transverse Chromatic
Aberration (LCA). LCA can introduce pronounced undesirable image artifacts as
shown later, partially since most lens systems correct for two spectra rather
than the three sampled by the sensor, due to complexity and cost limitations
[55]. Since the effect of LCA cannot be minimised in practice unlike with ACA,
this work will focus on addressing the intrinsic fidelity loss due to LCA,
exclusive of other types of aberration.

In order to fully differentiate both LCA and ACA types of chromatic
aberration, they are subsequently defined; other forms of aberration are
monochromatic and outside the scope of this work and are not considered
further. Chromatic aberration is the result of the refractive index of the
lens medium (typically some form of glass) varying with the transmitted light
wavelength, termed dispersion [69]; this is the definition of prismatic
behaviour.

With ACA, the focal point varies along the optical axis with the light’s
wavelength. As a consequence, the focal length of an object cannot be exactly
coincident in all three image planes, but only approaching this crucial ideal.
As light strikes the sensor plane, out of focus rays contribute to a circle of
confusion or bokeh, which manifests as a subtle coloured halo around the
boundary of an object in more extremal circumstances such as at the lens’s
widest aperture setting. The resulting introduction of image artifacts is
minimised as the lens aperture is stopped down or reduced due to the increase
in depth of field bringing the axially-misaligned focal points nearer. Many
modern digital cameras when in fully-automatic mode, balance the aperture size
preventing significant spatial frequency loss due to photon diffraction [66]
and a shallow depth of field, increasing focus selectivity and focus error as
a side-effect. Due to the camera’s automatic behaviour, ACA is nominally
minimised to imperceivable levels. Further, photographers generally reserve
the use of large apertures for shooting portraiture and low-light pictures
where maximum light flux is needed; lenses suitable for these situations with
larger maximum apertures, typically have stronger compensation for ACA.

Figure 1.1 shows the focal lengths of the primary wavelengths at different
displacements along the optical axis due to ACA. Walree [94] gives a clear
example of this, showing the effects of ACA when increasing the focal depth
sensitivity by means of opening the aperture, given in Figure 1.2. Due to the
ability to minimise ACA with lens parameters at picture capture time, this
typically does not pose a problem in practice.

For aforementioned reasons, ACA will not be considered further in this work.
This research will address the more nascent issue of LCA within digital
photography.

singlet elementoptical axisimage plane Figure 1.1: Axial Chromatic Aberration

Figure 1.2: Axial (left) and Lateral (right) Chromatic Aberration, adapted
from Walree [94]

###  1.3 Lateral Chromatic Aberration

LCA is the relative and non-linear displacement of the three colour planes,
across the image plane, conceptually shown in Figure 1.3. The result of this
chromatic distortion introduces undesired artifacts and thus information,
which leads to a perceived detail loss in the image, due to the misaligning of
coincident features throughout the image. Less perceptible impact in lower-
contrast areas reduces texture detail and generally tends to reduce the
perception that LCA compromises image quality.

LCA occurs along the sagittal or oblique direction to the optical axis. Figure
1.2, a top left image crop, illustrates that LCA only affects detail running
orthogonal to the line through the centre of the image, due to the different
distortion paths for the red, green and blue planes.

singlet elementoptical axisimage plane Figure 1.3: Lateral Chromatic
Aberration

Taking an example shown in Figure 1.4, shot in near-ideal real world
conditions with reduced aperture, no image shake and correct focus, we see
crops from various zones of the image showing significant LCA, estimated at
around an order of magnitude higher error than the Nyquist limits of the
sensor, i.e. up to around 10 pixels. The frame was taken at the optimal
aperture of f/8.0 at 18mm focal length (equivalent to 27mm on a full-frame
sensor) on a Nikon D90 digital SLR with a Nikon 18–200mm VR lens, considered
as basic professional equipment. It is clear that chromatic aberration limits
image detail, even without spatial frequencies approaching the Nyquist [84]
limits of the Bayer sensor array. Since the crops show the chromatic
aberration being asymmetric, LCA is therefore the dominant aberration and
attenuating quality. ACA occurs symmetrically to the image feature and is
shown to be far less significant.

Figure 1.4: Cropped sections exhibiting LCA

###  1.4 Central argument and research objective

The scope of this work is the correction of LCA using typical camera imaging
systems, and thus the most common lens designs and optical geometries, such as
large aperture, wide field-of-view angle and high zoom ratio goals; this can
be achieved only in digital post-processing steps. It is assumed that LCA can
be sufficiently addressed by the generalised distortion model. Correction for
optical decentering is discussed later. Though this would allow analysis and
correction with tilt-shift lenses, this is outside the scope of this work, as
these are comparatively infrequently used. Further, careful consideration will
be employed to avoid introducing additional user or system constraints, or
assumptions about the underlying system; this is a requirement for a robust
system.

Consequently, the primary objective is to develop a robust, generalised
algorithm able to perform across all real-world shooting conditions, and
without manual adjustment or input. Secondary objectives are to develop a
method to tangibly quantify image quality, and demonstrate the quality
improvement due to CA correction. Subsequent to this, analysis will be
conducted with synthetic images, and later with a typical imaging platform
with real-world images. Lastly, this implementation will be compared to
another widely-used solution of manual correction using Adobe Photoshop CS4.
Algorithms presenting a new correction technique and an algorithm to analyse
image quality improvement will be developed and evaluated throughout this
work.

###  1.5 Contributions

In this work, the minimisation of LCA by means of development of a robust and
automatic algorithm is assessed as a novel contribution in this field. An
algorithm to quantify distortion due to LCA is proposed and developed,
allowing determination of the optimal correction parameters and measurement of
recovery of fidelity loss. This is applied and evaluated to understand
algorithm robustness and correction performance. Methodology and workflow will
be demonstrated showing the ability to automatically correct real-world
pictures, allowing a closed-loop solution. Hence, these research findings are
expected to directly benefit those who are interested in this area at large.

In the next chapter, the context for the work will be presented with concepts
and definitions. Then, contributed work will be assessed relative to the
objective of this work.

###  1.6 Overview of methodology

The methodology that will be presented is characterised by the scheme shown in
Figure 1.5. Overall, there are two phases of operation: firstly, the
coefficient recovery phase, where correction data is computed from images.
Secondly, this correction data can be used to correct previously unseen
images, using data extracted from the database, resulting in images corrected
for LCA.

recovery phasecorrection phase source pictures coefficient recovery parameter
database correction corrected pictures source pictures Figure 1.5: Image
correction method overview

###  1.7 Thesis layout

Chapter 2, Background introduces concepts and equations to be used in later
chapters and looks at related constraints, presenting complete definitions.
Related work in this field is grouped and analysed relative to the objective
of this work — robust correction of LCA. Existing works are revealed according
to Chapter 3, Proposed method which presents high and low-level details of the
algorithms specialised to LCA correction, and breaks down each of the
processing steps. Chapter 4, Methodological evaluation defines technique to
show steps taken to validate increasing levels of capability. Finally, Chapter
5, Conclusion and further work discusses research findings, along with areas
for future consideration.

##  Chapter 2 Background

###  2.1 Introduction

This chapter aims to present major concepts and information relevant to the
understanding, measurement and correction of Lateral Chromatic Aberration.
These concepts will be used in later chapters.

The chapter is split into two important intertwined parts. To create better
understanding of the research topic, the first part will provide general
knowledge relevant to LCA, which includes optics, paraxial optics, lens
systems, optical properties, sensor properties, sensor blooming, image EXIF
data and parameters, lens variance, correction overview and test images. To
highlight the need of this research, the second part will review existing
literature on LCA and areas related to this. There are two types of image
correction which will be encountered: in-lens and post-processing. Post-
processing correction will later be examined in more detail, since the former
is only useful if control is available at the lens physical level. Post-
processing affords the freedom to arbitrarily change characteristics, required
for the intended use of this research. Subsequently, the knowledge gap found
will be discussed.

###  2.2 Optics

The field of optics encompasses a broad group of fields based around
transmission, modulation and sampling of electromagnetic spectra not limited
to visible light.

Concepts and equations directly relevant to this research are subsequently
introduced, and later employed.

####  2.2.1 Paraxial optics

Paraxial (lit. parallel-axial) optics model light as traversing optical
elements within an optical system with an angle of incidence from the optical
axis considered small [2]. Due to this low-angle approximation, wavefronts
propagating through the optical system are modelled as spherical chords, thus
radiating equally at all angles therefore free from aspherical aberration [2],
as illustrated in Figure 2.1.

aberratedwavefrontnon-aberratedwavefront Figure 2.1: Paraxial approximation

As the wavefront is modelled as a two-dimensional chord of a sphere, the
standard trigonometric $\sin$, $\cos$ and $\tan$ functions are used. The known
Taylor-series expansions for these are:

$\displaystyle\sin\theta$ $\displaystyle=$
$\displaystyle\theta-\frac{\theta^{3}}{3!}+\frac{\theta^{5}}{5!}-\frac{\theta^{7}}{7!}+\dots$
$\displaystyle\cos\theta$ $\displaystyle=$ $\displaystyle
1-\frac{\theta^{2}}{2!}+\frac{\theta^{4}}{4!}-\frac{\theta^{6}}{6!}+\dots$
$\displaystyle\tan\theta$ $\displaystyle=$
$\displaystyle\theta+\frac{\theta^{3}}{3}+\frac{\theta^{5}}{5}+\frac{\theta^{7}}{7}+\dots$
(2.1)

Discarding orders higher than one in Equations 2.1, the approximation becomes:

$\displaystyle\sin\theta$ $\displaystyle\approx$ $\displaystyle\theta$
$\displaystyle\tan\theta$ $\displaystyle\approx$ $\displaystyle\theta$
$\displaystyle\cos\theta$ $\displaystyle\approx$ $\displaystyle 1$ (2.2)

Applying this to the lens model, $\theta$ is the ray angle from the optical
axis; lens systems with these assumptions are termed first order optics, since
they are valid up to only the first order approximation.

####  2.2.2 Lens systems

A lens is defined as one of more groups of lens elements through which light
is focussed. Single Lens Reflex (SLR) camera lens systems typically have four
to twenty lens elements in a number of moveable groups [49]. Lens elements are
grouped to minimise dispersion, which occurs where different wavelengths of
light are refracted through the medium at angles varying with distance from
the optical centre. Since dispersion is a physical property of the material,
it cannot be eliminated. Lens designers strategically employ lens elements
with multiples types of glass, with intentionally differing Abbe numbers [55],
also known as dispersion or refractive index. When elements are placed in a
pair, cemented back-to-back, the number of medium transitions is reduced and
thus light loss due to reflection; significant dispersion angles can be
avoided. Such an arrangement is termed a doublet and illustrated in Figure
2.2. The second element in the pair is often selected to disperse with the
opposite sign of the first. Flint and crown glass types are most commonly used
[50] for this purpose. Techniques of applying an anti-reflective coating are
used on selected elements to minimise additional light loss at medium
transitions [73].

crownflintelementelement Figure 2.2: Dual thin lens (“doublet”) with CA
correction

From these implementation constraints, three resulting compromises occur:
firstly, there are material choice constraints, secondly, at worst there are
twice as many medium transitions for propagating rays with using doublets, and
lastly, there are additional manufacturing requirements such as precision and
surface geometry, where the pair of elements are interfaced together forming a
doublet. An additional set of issues arise as a result: there is increased
light loss from the more frequent medium transitions, secondly, there is also
more light scattering from surface reflection [39], and a single optimal
material cannot be selected alone, due to uncontrolled dispersion.

In order to understand light loss due to surface reflection, the angle of
refraction can be computed using Snell’s law from Figure 2.3 and:

$\theta_{i}$$\theta_{t}$$n_{i}$$n_{t}$ Figure 2.3: Snell’s law

$\frac{n_{i}}{n_{t}}=\frac{\sin\theta_{i}}{\sin\theta_{t}}$ (2.3)

where $\theta_{i}$ is the angle from the surface normal of the incident ray,
$\theta_{t}$ is the angle of the traversing ray, $n_{i}$ is the refractive
index of the primary medium, and $n_{t}$ is the refractive index of the
secondary medium.

Using fluorite lens elements was later found to minimise CA due to its lower
refractive index than other types of glass [73]. Since fluorite is brittle and
requires significantly more mechanical work to grind, it features in typically
one lens element in high-end lenses due to the significant cost premium.
‘Ultra-low Dispersion’ (UD) lens elements were developed as a way to get lower
refraction by chemically doping normal glass to lower the refractive index, at
a lower cost; this allowed entry into semi-professional lenses, where one or
two UD lens elements are incorporated to reduce CA. In more recent times,
special lens elements using diffraction grating on stock glass have been
employed [92]; however, they require premium and expensive manufacturing
techniques only more recently available [46]. The cost of lenses incorporating
such developed elements is beyond consumer reach at large.

Today, the use of diverging-converging doublets termed achromatic doublets,
shown in Figure 2.2, remains the most ubiquitous way of minimising CA [86],
yet this is known to be an imperfect solution due to three reasons: firstly,
the doublet optimises correction for two primary wavelengths of red and blue,
around the green wavelength, leaving an uncorrected secondary spectrum [50].
Secondly, cheaper and/or simpler lenses are modelled using paraxial and thin
lens approximations, termed Gaussian optics, which are suitable for
applications where quality is not the primary concern [50]. These modelling
steps result in first-order residual error. More sophisticated lens designs
employ computational ray-tracing through the lens system, such as offered in
some commercial lens design software, such as Sinclair Optics’ OSLO [71].
Lastly, the number of interfacial surfaces is therefore greater (four, rather
than three for a doublet) as previously mentioned, increasing light loss.

####  2.2.3 Lens parameters

With the construction of lenses, a number of element groups are designed to
move with manual linkage or electrical control. Two parameters are observed
from element movement, and a third invariant of element movement. Firstly, the
lens focal length relates to the magnification of the image. Secondly, the
focussing distance allows selection of the focus plane. Lastly, the aperture
size can be changed to vary the depth of field. All three adjustments affect
the path of light, so are recorded in the image metadata, thus available.

####  2.2.4 Optical properties

Due to the simplicity and consequent limitations of paraxial optics, a more
complete framework to model lens systems was sought. In 1956 Ludwig von Seidel
developed a model defining various types of orthogonal aberrations [95],
allowing aberrations to be decomposed into constituent types and expressed
independently.

The complete Seidel equation [54] using Zernike wave-front notation [97] is:

$\displaystyle W(r,\phi)$ $\displaystyle=$ $\displaystyle w_{020}\ r^{2}$
$\displaystyle+$ $\displaystyle w_{040}\ r^{4}$ $\displaystyle+$
$\displaystyle w_{131}\ \eta\ r^{3}\cos\phi$ $\displaystyle+$ $\displaystyle
w_{222}\ \eta^{2}\ r^{2}\cos^{2}\phi$ $\displaystyle+$ $\displaystyle w_{220}\
\eta^{2}\ r^{2}$ $\displaystyle+$ $\displaystyle w_{311}\ \eta^{3}\ r\cos\phi$

where $\eta$ is object size, $r$ is displacement from the optical centre,
$\phi$ is the angle clockwise from vertical of the ray, and $w_{ijk}$ is the
coefficient describing the magnitude of aberration. Where only distortion
aberration is present, a number of terms become zero, leaving:

$W=w_{311}\ \eta^{3}\ r\cos\phi$ (2.5)

Following on, it is known that the cosine function can be written as an
infinite Taylor series:

$\cos\phi=1-\frac{\phi^{2}}{2!}+\frac{\phi^{4}}{4!}-\frac{\phi^{6}}{6!}+\frac{\phi^{8}}{8!}-\dots\\\
$ (2.6)

This becomes intractable, so a finite number of terms in the Taylor series are
used, approximating the cosine function, here up to the fourth order:

$\cos\phi\approx\phi-\frac{\phi^{2}}{2!}+\frac{\phi^{4}}{4!}\\\ $ (2.7)

The Seidel approximation proves the ability to mathematically model a set of
lens aberrations with a Taylor series. Since this work applies to general,
aberrated lenses, it is not possible to make assumptions about lack of
particular aberrations, thus the only solution is to use the generalised
Taylor series:

$\displaystyle r_{dest}$ $\displaystyle=$ $\displaystyle
a{r_{src}}+b{r_{src}}^{2}+c{r_{src}}^{3}+d{r_{src}}^{4}+\dots$

where $r_{src}$ is the image sensor plane, and $r_{dest}$ is the corrected
image plane. $a$ is the overall scaling factor, and, $b$ and $c$ exert
influence on the non-linear chromatic distortion of the image plane. This can
be used to non-linearly register one plane against another.

####  2.2.5 Sensor properties

In order to develop an imaging mechanism and thus a sensor that minimises
information redundancy, Bayer [5] developed the concept of a Colour Filter
Array (CFA) optimised to the spectral sensitivity of the Human Visual System
(HVS). It is known that the eye has three photoreceptor cones with peak
sensitivity at three wavelengths, 564–580nm long, 535–545nm medium and
420–440nm short, and twenty times the number of luminosity-sensitive receptors
or rods [102], sensitive to a wider spectral range and thus suited to low-
light level vision, as they absorb wavelengths across a far greater band.
Earlier research conducted by the Commission Internationale de l’Eclairage
(CIE) showed that the eye is more sensitive to the medium wavelength than the
low or high, leading to the development of the photopic weighting function
$V(\lambda)$ [18, 48]:

$\displaystyle Y$ $\displaystyle=$ $\displaystyle V(\lambda)$ (2.9)
$\displaystyle=$ $\displaystyle 0.299R+0.587G+0.114B$

where $Y$ represents luminance of a pixel, $V$ is the photopic equation,
giving the HVS response to light of frequency $\lambda$; $R$, $G$ and $B$ are
the pixel values in the red, green and blue planes respectively. The Bayer
array sensor thus maintains a quadrant grouping of elements, where two are
sensitive to the medium wavelength, one to the short wavelength and one to the
long wavelength. Specialising the elements to a particular light wavelength is
achieved by a set of colour filters directly in front of the sensor. One
negative side effect of filtering the light is that light flux illuminating
the sensor elements is reduced by at minimum 66% [51].

Processing steps are needed to reconstruct the CFA to the typical arrangement
of stored image data of red, green and blue image values per pixel. Simple
linear interpolation, shown in Figure 2.4, is known to leave chromatic edge
artifacts [4], so further algorithms have been developed such as Variable
Number of Gradients (VNG) [15], Patterned Pixel Grouping (PPG) [44] and
Adaptive Homogeneity-Directed (AHD) [42] interpolation.

imaging sensorimage planesinterpolation Figure 2.4: Bayer demosaicing
principle

Other sensor types exist, such as the Foveon X3, which uses a stacked array of
photodiodes [33]. Since the CFA is ubiquitous, other sensor types are not
considered in this work. The algorithms in this work equally apply and
theoretically benefit from the higher chromatic spatial frequency information,
allowing more accurate LCA recovery.

####  2.2.6 Sensor blooming

The term blooming applies exclusively to Charge Coupled Devices (CCDs) due to
their implementation. This describes a specific situation that occurs: the
charge caused by the potential stored at a CCD element overcomes the
dielectric between one of more surrounding elements, and charge is transferred
away or leaked. The manifestation of this is typically that purple halos
appear around areas which have experienced excessive saturation due to there
being twice as many elements sensitive to green per unit area as blue or red.
The mechanism is depicted in Figure 2.5.

imaging sensorover-saturatedelements Figure 2.5: CCD blooming principle

Other types of imaging sensor, such as Complimentary Metal Oxide Semiconductor
(CMOS) sensors are implemented in a different way and indirectly avoid this
mechanism [30].

###  2.3 Image EXIF data and parameters

Image meta-parameters are stored in the EXIF [26] section in the image source
file. This contains parameters which are not needed to correct the image but
are essential to reusing the correction data later and secondarily storing
useful state once correction is performed.

The camera body identifier and lens identifiers are used to select the right
set of correction parameters; the focal length, aperture value and focus
distance are affected by the position of the lens element groups, thus need to
be known to find the nearest neighbour set of correction coefficients. The
camera body orientation and image size information is needed to later store
the lens decentering values, should this be calculated.

Subsequent to LCA correction analysis, a LCA corrected flag can be set to
indicate not to perform LCA correction on the image again, and the correction
coefficients can be stored in the image file and database for later use.

[botcap,caption=Example stored EXIF tags and image parameters,label=image-
params]ll From camera Lens identifier Nikon AF-S DX VR 18–200mm Camera body
identifier Nikon D90 SLR Focal length 18.0mm Subject (focus) distance 3100mm
Lens aperture value f/8.0 Image orientation Top-left Image width 4352 pixels
Image height 2868 pixels Calculated Computed image X-centre 2172 pixels
Computed image Y-centre 1432 pixels LCA corrected False LCA R-G correction
coeff 1 0.99855155 LCA R-G correction coeff 2 0.0032236 LCA R-G correction
coeff 3 -0.00190334 LCA R-G correction coeff 4 0.00120401 LCA B-G correction
coeff 1 1.00113412 LCA B-G correction coeff 2 3.02337291e-4 LCA B-G correction
coeff 3 -2.45039357e-3 LCA B-G correction coeff 4 1.76123214e-3

###  2.4 Lens variance

Due to variance in manufacturing processes, there are sub-micron imperfections
in lens element seating and alignment [87]. A tight Gaussian curve of
variation is seen with quality control tests causing rework of lenses outside
the tolerances. Once in consumer hands, wear in the guide grooves for the lens
mechanics, particularly for zoom lenses, increases sensor data error.
Ultimately, this causes lenses to slowly diverge from the original
specification. Professional users have reported variation among particular
lenses, with some units attenuating spatial frequency more than expected [17].

This field variance presents a key requirement: the correction data from one
lens unit may not optimally correct another lens unit of the same model, so an
adaptive technique which optimises for the behaviour of the lens from which
images were taken is needed.

###  2.5 Correction overview

Correction of LCA with post-processing steps in a digital workflow is, at
present, the only solution offering complete correction. Therefore, post-
processing steps are slowly being incorporated to indirect (in-camera)
processing [22], since it is simply not possible to design out LCA in a lens
with more than one element. Various techniques have been developed from manual
adjustment of dual-linear correction parameters using modern image processing
software [7], to more elaborate and calculated step-wise procedures [57]. Both
approaches are error-prone due to the user’s subjective observation feeding
back into the adjustment process, resulting in a tendency to be inaccurate and
slow. Moreover and critically, the work to correct a single image does not
contribute useful correction values towards an accurate model of the lens due
to lack of information re-use. Indeed, manual LCA correction by inspection
relies on the user assessing the overall resulting picture and is inherently
based on psychovisual feedback, making this an unsuitable process for accurate
and consistent data capture. Areas which are less perceptually apparent, e.g.
lower spatial-frequency textures, or lack of contrast edges may be under-
corrected.

The second key limitation post-processing introduces is that of complexity. A
popular correction mechanism has two variables to manually adjust, which
control linear scaling of the red plane relative to the green plane, and the
blue plane relative to the green plane. Due to non-linearity in the lens
system, this linear mapping is an approximation. Even with more complicated
methods, accuracy is still limited and should be at the sub-pixel level.
Additional to this, if LCA correction is performed before demosaicing the
sensor array, the advanced demosaicing algorithm will be able to use the true
edges and not edges introduced from LCA artifacts. Further, advanced
demosaicing algorithms can introduce false colour aliasing [56], which would
change the true LCA-free minima. Non-linear interpolation, sub-sampling or
sharpening may change spatial information in an undesirable way, such as
emphasising high-frequency noise. To prevent this, the image should be
converted from gamma corrected to linear before interpolation and back after
demosaicing.

A number of camera manufacturers (e.g. Nikon [47]) have recently introduced
software LCA correction internal to some of their medium-level digital SLR
cameras. This provides a practical integrated solution to tackle the issue of
significant LCA on wide-angle lenses, and is a value-adding feature touted to
improve image quality, giving an edge over competing products. Nevertheless,
there are compromises which restrict the usefulness of this in practice.
Firstly, when capturing maximum sensor detail for off-line processing,
photographers shoot in raw image format, so in-camera LCA reduction is
therefore not performed; it is only available when capturing in lossy JPEG
format. Additionally, since processing is implemented on a minimal embedded
platform, tradeoffs and approximations are employed to give reasonable
processing time, and will be weighed against other tasks performed at capture-
time, including picture analysis, white-balance correction, Bayer
interpolation and JPEG compression. Most implementations use some form of
Discrete Cosine Transform fast-integer approximation or other techniques to
speed up compression encoding [104]. Lastly, quantisation tables are used to
truncate encode detail to the desired level.

###  2.6 Test images

The input images used in this work are selected intentionally to present a
number of typical shooting conditions. The input images are shown later in
Figure 2.6. Images 1–8 are shot outside on an overcast low-brightness day,
where a small aperture around f/8.0 would be used to minimise ACA. ISO
sensitivity is higher at 400, and exposure time is longer at $5\times 10^{-2}$
to $1\times 10^{-2}$ seconds. Images 9–12 are typical light-constrained indoor
shots with the aperture balanced around f/6.3; an increased ISO around 800 is
used and much longer exposure time in the range $5\times 10^{-1}$ to $1$
second. Subsequent images 13–26 are shot in clear and bright conditions,
allowing unrestricted use of shooting parameters, thus using an f/11.0
aperture for depth of field, higher corner detail and less ACA; ISO is set to
200 and exposure time between $1\times 10^{-3}$ to $2\times 10^{-2}$ seconds.

Additional to the real-world images, a CA-free picture is rendered using
POVray with a high-complexity office scene [74], shown later in Figure 4.3.
Further, simple chequerboard patterns are rendered for basic validation, shown
later in Figure 4.1.

From visual inspection across a number of the test images, the most visible
LCA artifacts (i.e. at high-contrast edges) are seen to be around 3–4 pixels
in displacement. Since the red and blue planes were linearly interpolated from
a Bayer array (due to simple demosaicing), there is around one pixel of
uncertainty.

Figure 2.6: Real-world test images

###  2.7 Literature review

Available research relating to LCA correction covers a wide spectrum of themes
across varying application environments from scientific measurement to
practical areas applicable to photography. Directly and indirectly relevant
works are identified and grouped thematically into optical correction,
typically achieved through lens design enhancements, and computational
correction (typically achieved through software post-processing) areas. The
contribution of the literature is then assessed relative to the goal of this
research. For completeness, CA is briefly addressed in the wider field of non-
camera optics relating to other imaging systems, including in the fields of
microscopy, diffractive systems and contact lenses. Finally, the knowledge gap
is identified and discussed.

####  2.7.1 CA and non-lens systems

Since CA affects waves of differing wavelengths propagating through media in
general, the effects of this are found in many non-lens systems (though in
principle, it affects all such systems). Much of the research in these areas
shows CA as having profound effects in three major fields outside the
classical use of optics: ocular CA affecting the HVS, electron microscopy and
diffractive systems.

Firstly, the study of LCA in the human eye is a buoyant field, driving
advances in laser surgery and other treatments. Hay et al. shows the HVS to
have intrinsic adaptation to CA [41], though LCA alone is shown to contribute
significantly in vision acuity loss in humans despite the fovea covering only
an inner 3 degree arc [88, 89]. CA introduced from contact lenses or
eyeglasses has been shown of higher significance, and a wide body of
literature is available on this. Powell introduces a contact lens design that
corrects ACA without increasing LCA [76]. Likewise, Zhang et al. finds the use
of dual lens elements has proven to address LCA without the previous increases
of ACA [105].

Secondly, correction for chromatic aberration in scanning electron microscopes
involves third order and greater wavefront correction, based largely on the
emergent work by Scherzer [80]. Later research by Tucker et al. used cubic
phase plates for wave-front coding to increase the depth of field in
microscopy while addressing ACA [91].

Lastly, the ability to use Fresnel and diffraction gratings with optics to
compute Fourier transforms is limited in polychromatic light sources due to
frequency-dependent variation; Ferriere and Goedgebuer developed a lens
configuration which eliminated chromatic aberration by using an exactly
symmetrical selection and placement of lens elements along the lens [29].

####  2.7.2 Lens-system optical correction

A group of the literature addresses correction or compensation at some level
of chromatic aberration by means of modifying the optical path. Since this
occurs physically prior to image capture, it requires in-loop design knowledge
and modification typically via hardware or design changes. Due to the less
flexible nature and the expense of implementing changes in this way, this area
lacks the accessibility that alternative correction mechanisms offer.

Grey [36] conducted analysis on the inter-dependency of the various types of
aberration within aberration theory, finding that adjustments to correct for
one type of aberration affected other types, due to the non-orthogonal nature
of optical aberrations. This work presented early conceptual methodology to
computationally refine lens design using error minimisation.

The earliest literature to directly address CA introduced the minimisation of
CA via active optical correction by Willson and Shafer [99]; this was achieved
by using three independent optical paths for the three primary wavelengths
through zoom lens motor micro-adjustment. This allowed an overall reduction in
LCA for a specific focal length. Since three exposures are needed, temporal
coherence is lost; this presents real-world impracticality with the motorised
zoom lens constraint. Due to three discrete image samples, this method is thus
unsuitable for scenes with movement.

In optimising lens design for lower CA (and other aberrations), Fang
introduced the use of genetic algorithms [28]. This addressed the complexity
associated with lens element selection from the wide variety of lens elements
available, and from a range of sources. He demonstrated a measurable reduction
in LCA and other optical aberrations.

Spatial Light Modulators (SLMs) have been developed in certain application
fields where cost is not sensitive, e.g. microscopy. Millan et al showed
application of this in diffractive optics and modifications to the light exit
aperture to reduce LCA [3].

More recent research has emerged showing ray-tracing through optical lenses in
forward and backward directions using lasers and CCD sensors to obtain
empirical chromatic aberration data [81]. This was shown to address the
complexity and assumptions in software lens models though iteratively
optimising the design of the lens system.

Indirectly related at the sensor capture level, CCD blooming presents
significant issues and is often misinterpreted as CA due to the chromatic
changes it causes. This is crucial to particular fields outside general
photography, such as in astrophotography: the contrast of planets illuminated
by the sun against free space presents very high contrast achievable in
adjacent CCD cells, saturating the charge buckets and causing charge leakage
to adjacent cells. At the CCD cell level, the introduction of additional
circuitry to mitigate this has been shown to reduce the impact of this issue,
at the expense of increased noise and complexity [70]. Notably, CMOS systems
are not fallible to this issue, so they remain a superior choice for high-
contrast optical systems.

In a wider field outside correction of CA, active optics have been in use for
some time. One of the most significant and earliest fields was to use the sub-
field of adaptive optics in connection with astrophotography. This is employed
to adapt the surface of a telescope mirror to approximately rectify the
incoming wavefront as far as possible. To achieve this, Hartmann-Shack
wavefront sensors and a tight control loop to update an array of actuators are
used [82]. Related to this technology, Digital Micromirror Devices (DMDs)
allow high-speed deflection (and therefore switching) of the optical light
path [43] within discrete pixels. DMDs have been widely deployed in consumer
image projectors.

Although not directly relevant to CA reduction, Willson addressed the
calculation of lens intrinsic parameters in motorised zoom lenses over a wide
range of lens parameters, identifying additional problems [100]. Later, he
went on to model the relationship between various extrinsic camera parameters
for a motorised zoom-lens system [98], based around Tsai’s model [90]. This
work proved that interpolation is successful in obtaining values at
intermediate lens parameters from a discrete set of parameter. Chen et al.
went further to show bilinear interpolation on lens parameters holds for
manual (i.e. non-motorised) lenses [16]. This allows measurements to be taken
for a number of discreet lens parameters, interpolating between available
parameters to correct for a tuple of parameters previously not seen.

####  2.7.3 Computational correction

While design changes, introduction or modification of hardware present cost
and engineering constraints, the modification of image data post-capture is
attractive due to the accessibility and freedom it brings. Correcting images
computationally is inherently in the post-processing phase, as the image data
is digitised after being affected by the analogue optics.

Relative to the optical correction group of works which can be considered
comparatively mature, this area is rapidly developing. Some methods within
software correction, including forward and inverse Fast Fourier Transforms and
deconvolution, introduce significant error at the typical quantum depths of 32
or 64 bits per pixel. Future image processing is thus reliant on the
introduction of GPU and processor architecture support for 128 and 256 bits
per pixel quantum depths to minimise introduced noise, which will allow far
more creative license, such as refocusing and near-total aberration correction
via high-precision iterative deconvolution.

Overall, the majority of software techniques found present quantification,
analysis and algorithms centred around fixed lens systems, i.e. for a tuple of
given lens settings. This is traditionally sufficient for most applications,
as for example web cameras commonly used in vision systems have a fixed focal
plane and aperture, however with increasing imaging requirements and
applications partially driven by continuously increasing specifications and
costs, this area has steeper requirements and expectations.

Early work in the field suggested the use of image warping to correct LCA [35,
79]. Boult understood solving this issue numerically had advantages and
introduced LCA correction via digital imaging with cubic splines as a function
of radial distortion [101]; this resulted in a similar error reduction as with
Willson’s active optics. An achromatic chess board target was used to measure
improvement.

Dersch started a project in 1998 to develop a toolkit to allow manipulation of
images for the correction of distortion aberration [25]; this relates to LCA
correction through the use of warping for aberration correction. This was
later extended for the purpose of automating image identification and
registration for panoramic photography and used Lowe’s SIFT algorithm [61] to
extract stable image features [11]. A similar and early approach in this work
used the SIFT algorithm to identify image feature points on each plane for
converging via non-linear registration, however instability and inaccuracy
were found due to lack of information from reliance on discrete image
features, which also limited radial coverage and resulting in under-correction
at image extremities.

Later, Dersch’s project was forked and openly developed. Watters identified it
as being useful for the correction of LCA, presenting a manual technique for
finding approximate LCA correction coefficients [96]. The steps consisted of:
the red and blue channels were manually aligned to minimise visible LCA at an
initially chosen radius, recording the pixel displacement and radius into a
spreadsheet. Next, this was repeated for four further radius intervals.
Lastly, the number of data points was selected, and the curve fitting function
in Microsoft Excel was used to generate correction coefficients. This
technique was subsequently refined by Krause to simplify the linear scaling of
channels by the development of a software tool [57], though the red and blue
channel correction coefficients were still curve-fitted via spreadsheet
software as before.

Notably, Park et al. used feature points for measurement and subsequent
elimination of barrel distortion on zoomable lens systems [72]. Benhimane
developed a similar approach to correct for distortion aberration via corner
detection against specific calibration patterns [6].

Remondino concluded that LCA was the dominant form of CA and thus needed a
robust correction mechanism [78], contextualising the basis for this work.

Kaufmann used EOS Systems’ PhotoModeler software to obtain a simple
translation and linear scaling for the red and blue channels relative to green
[53], citing the translation was necessary to compensate for slight sensor
and/or lens element misalignment. Difficulty was reported in detection of
false-positive features when calibration was attempted in real-world scenes,
concluding in the suggested use of bespoke achromatic calibration targets.

The work of Heuvel and Verwaaal showed that application of available camera
barrel distortion correction software—in this case EOS Systems’ PhotoModeler
on a per-plane basis—was largely unsuccessful with their fixed fish-eye lens
application [93]. Least squares minimisation was performed to reduce
registration standard deviation of the 90 targets they measured down to 0.15
pixels. This left large registration error cited as up to 6.0 pixels, which
was attributed to other issues, such as CCD blooming despite using linear
plane scaling approximation. Further tests were conducted using manual linear
scaling to reduce image artifacts observed, using Picture Window Pro 4.0,
citing “significant visual improvement”, but it still proved an incomplete
solution [93].

A self-calibration model proposed by Cronk used coloured circular features to
calculate intrinsic parameters [23]. Rudimentary blob-detection was used with
five differently coloured features and multiple views. Quantification of the
calibration results showed sub-pixel registration across much of the image,
however accuracy was proven problematic citing a worst-case error of “several
pixels”.

Mallon et al. introduced fully automated correction of LCA through post-
process warping the three colour planes from images of classic monochromatic
chess-board calibration targets [63]. Iterative minimisation was achieved
through analysis of the colour histogram from the calibration target.

It has been found that variation among zoom lenses due to wear and hysteresis
is ultimately significant for accurate correction mechanisms [17], challenging
these earlier models.

Luhmann [62] used the method of self-calibrating bundle adjustment on features
from a calibration pattern, using interior and exterior parameter correlation;
correspondence was introduced between channels correcting for LCA. This was
achieved with per-channel feature measurement with opaque software Ax Ori
Axios 3D. Fryer and Brown [32] extended this to use plumbline calibration for
calibration without exterior orientation. An improvement to feature-point
accuracy of factor 1.6 was shown. Later, he demonstrated the effect of LCA via
channel subtraction [62], showing that it can be measured in a continuous way
to produce an edge-map.

Remondino [78] et al. proposed a method employing the use of colour filters to
selectively calibrate individual colour planes with an achromatic calibration
target and suggested the use of capturing the raw sensor data. They were able
to show the significance of LCA with a maximum registration error of 10 pixels
with a second-order radial distortion model based around earlier work [90].
This was shown for four different focal lengths and varied with changes in
focal length.

Taehee used a technique of edge refinement to minimise LCA [60]; a filter was
applied in both directions to adjust pixel colour when meeting various colour
saturation conditions. The behaviour of this approach was not analysed in
pictures where LCA may be falsely be measured, leading to the conclusion that
it may not be robust. Moreover, this technique does not account for
geometrical changes introduced by LCA, and adapts local regions.

The recent work of Cecchetto [14] used inter-channel feature point
correspondence, with statistical outlier filtering to correct LCA. As
previously found, this technique highlights non-robust behaviour due to
relatively lower-density of discrete feature points.

Pomaska superficially showed that subtracting per-channel pixel values in
Photoshop could be used to visualise changes in distortion (albeit for barrel
distortion) [75]. This proves the mechanism of allowing a continuous approach,
as opposed to using discrete image features for manual correction of barrel
distortion. This allows the key technique of continuous correction to maximise
information utilisation.

###  2.8 Existing approaches

The overriding advantage of the existing approaches grouped into the design
and manufacture improvements is that the solution becomes intrinsic to the
lens. Therefore, correction cannot be omitted later, and there is no potential
for spatial frequency loss due to aliasing. Further, once the changes are
incorporated into the lens design, there is no further time or production
costs, and material costs from selection of more exotic glass materials with
lower dispersion are known up-front.

In research grouped into image post-processing, the primary advantage is
compelling: end-users or consumers of the optical hardware would be able to
perform this or similar calibration and/or correction in order to minimise the
effect of LCA and using arbitrary hardware. This has two issues though,
firstly digital correction has a steep requirement of potentially a complex
number of steps, and secondly that developed computational models and
algorithms are needed. Optical or design techniques intrinsically avoid these.

###  2.9 Summary

Evaluating literature grouped into computational correction, a number of
research works perform prior calibration against either the classic chess
board or a target bespoke to the application. This introduces challenges in
designing and fabricating a suitable target, error introduction from target
mishandling (e.g. a printed target being imperfectly planar), and the need to
perform calibration for each individual lens and parameter set.

The intersection of LCA correction and camera optics shows a narrow field,
mainly understood due to LCA being an increasing concern with the availability
of high-resolution sensors and reasonable fidelity optics. The application of
a number techniques in other imaging areas has been tested against the
application of LCA correction and highlights some areas of non-robust
behaviour. This is particularly the case where post-capture techniques for LCA
correction are manual and based on visual judgement; this demonstrates a wide
knowledge gap with automatic or closed-loop correction.

Consequently there is a clear opportunity to develop techniques tailored for
the micro and interrelated adjustments needed for stable LCA detection and
correction to a sub-pixel level. A model that takes advantage of the full
continuity of all the pixels on the image is proposed, rather than using a
comparatively small subset of the information available from discrete feature
points.

The next chapter will deal with recovering longitudinal chromatic aberration
distortion information from a collection of images, via a completely automated
workflow. Using appropriate image decoding and analysis, accurate chromatic
aberration correction can then be applied, exceeding the resulting quality of
most hardware and solutions and with no requirement for user-expertise and
allowing totally automatic correction.

##  Chapter 3 Proposed method

###  3.1 Introduction

In this work, an algorithm that eliminates LCA is proposed. This is achieved
by compensating for the prism-like deflection of light through the lens
elements used in capturing a photograph. Since the camera’s sensor filters
light down to three discrete wavelengths rather than a continuum of
wavelengths, it is proposed that the two sets of image data (planes)
reconstructed from red and blue elements of the Bayer CFA, will undergo non-
linear radial warping to map them onto the third (green) plane.

As the camera’s sensor measures the luminance of light at these wavelengths
across the same plane, separation of the sensor cells with the same
wavelengths into different planes gives three incomplete arrays; typically
cells are spatially divided into groups of four as shown earlier in Figure
2.4: two with a green cell-filter and the other two with red and blue filters.
This provides the most cost and complexity-efficient implementation due to the
eye’s peak sensitivity near the 555nm (green appearance) wavelength [27],
thereby giving increased spatial resolution in the green plane. The CIE 1924
colour appearance model [18] uses this basis. To derive an image plane with
the quarter or half per-colour pixel occupancy from the sensor, linear
interpolation is used and results in a continuous plane. Since the green plane
has twice as many samples per unit area than red or blue, the red and blue
planes are selected for later distortion to register against the green plane;
this minimises error induced from pixel aliasing which would reduce the
maximum spatial frequency of the resulting image.

###  3.2 Approach

The proposed algorithm remaps pixels according to a radial distortion model
from the source plane; this is done to align luminance features and detail
with that of the pixels in the green plane, removing LCA. This is feasible
through exploiting inter-planar luminance correspondence, for example due to
objects commonly contributing luminance information over multiple RGB planes;
even green leaves on a tree have some luminance detail in the red and green
planes, but moreover their background will be visible and give strong edge
detail. Geometric differences are thus observable [62], allowing an error
function to be derived and subsequently minimised to give inter-plane
congruency.

From Snell’s law, differing wavelengths will refract at different angles
through a medium, thus the resulting geometrical error is radial from the
optical centre. As shown earlier, modelling this radial distortion using a
Taylor series as shown in Equation LABEL:eqn:seidel-approx and up to the
fourth order is suitable. The resulting polynomial allows a remapping function
to be generated, used to compute radial pixel displacement from the optical
centre; this is used to reverse the chromatic distortion due to LCA.

Figure 3.1: Plane subtraction first observed by Luhmann

Since the internal lens geometry and precise positioning of lens elements and
aperture are abstracted somewhat by imprecise sensor measurement and thus
unknown, the function of how light will distort cannot be known or computed a
priori; manufacturers typically ray-trace through theoretical lens models, and
equations used do not apply to end-use of the system [58]. Lack of insight and
complete specification of internal lens geometry and construction prevents
direct mathematical or computational methods of finding the radial distortion
model coefficients providing least error. Instead, an error function is
defined and minimised iteratively. A novel technique is employed for the error
function: for each R-G and B-G pair, iterate over all the pixels in the image
planes, summing the absolute difference between each plane; a similar
technique was first observed by Luhmann, though using subtraction [62], shown
in Figure 3.1. As LCA introduces spatial incongruities between an object’s
pixels in different planes, the absolute difference is numerically lower when
there is complete alignment. Typical visible incongruities are seen in images
as colour fringes around high-contrast features, such as dark tree branches
against the light sky, and particularly towards the image extremities where
the aberration error is typically maximal, and zero at the optical centre from
Seidel’s aberration theory, seen in Equation LABEL:eqn:seidel-equation. An
appropriate non-linear error minimisation function is employed to adjust the
Taylor series coefficients; the distortion is applied against the image plane,
and the error function executed and further iterations performed until
sufficient accuracy is met. This technique is only valid in the untransformed
RGB colourspace; conversion to other colourspaces samples values from multiple
channels at coincident points, defeating the separation needed to preserve the
discrete chromatic information. The application of this method is shown in
Figure 3.2.

Figure 3.2: Cropped sections from red-green uncorrected absolute difference
map

Since the optical path through the lens system varies with focussing distance,
focal length and aperture size, the corresponding tuple of correction
coefficients as applied to the polynomial, cannot be applied to pictures with
non-identical lens parameters. For practical use, a database of the lens
parameters is constructed allowing applicable correction coefficients to be
found for any given set of lens parameters; this is chosen since correction
parameters for a given tuple of lens parameters can be held as state for rapid
later use. It is hypothesised that a significant number of pictures across
numerous lens settings need to be optimised and correction parameters added to
the lens database, to give sufficient data for later high-fidelity offline
correction on unseen images.

Since the computational cost of recovering the LCA correction coefficients is
high, and a single remapping correction pass is comparatively low-cost, the
overall procedure is performed in two stages: coefficient recovery and image
correction, shown previously in Figure 1.5. This allows efficient re-use of
the correction coefficients.

####  3.2.1 Error function

The function which is minimised to undo the distortion that LCA has introduced
is core to the algorithm. It operates across all pixels in the source plane
buffer and distorted plane buffer, returning a scalar value; lower values
indicate less LCA incongruity. The function is:

$\displaystyle e(P,M)=\frac{\sum_{P\in M}\mid P_{g}-P_{w}\mid}{n}$ (3.1)

where $P$ is a pixel in the set of pixels in mask $M$, of which $P_{g}$ is
pixels value in the green plane and $P_{w}$ is the value in the working plane
(i.e. the red or blue plane) and $n$ is cardinality of set $P$.

####  3.2.2 Decentering calibration

Lens decentering relative to the image centre is useful for the benefit of
completeness according to Brown’s method [10] and is not used as part of the
algorithm presented; it is simply integrated by changing the distortion centre
offsets, thus it is intrinsic to the methodology.

The offsets due to optical decentering are relatively small with respect to
the magnitude of LCA aberrations seen with production SLR lens hardware,
itself by definition non-existent at the optical centre. Any decentering
offset outside the central image zone would have a diminishing effect as the
distance from the optical centre increases [52], therefore the impact optical
decentering has on image quality is negligible and is not considered further.
Existing works using manual linear adjustment omit this, though it could be
considered in a future step.

###  3.3 Coefficient recovery

1

2

3

input : Image to be corrected complete with metadata

input : Database of lens intrinsic parameters

output : Database entry for correction coefficients at given lens parameters

output : Optionally, corrected image with updated metadata

4

5$\textnormal{{Image}}\leftarrow$ LoadImage(__) ;

6 $\textnormal{{Metadata}}\leftarrow$ Image.Metadata ;

7
$\textnormal{{Intrinsics}}\leftarrow\textnormal{{Database}}[\textnormal{{Metadata}}.\textnormal{{Lens}}].\textnormal{{Intrinsics}}$;

8 $\textnormal{{Planes}}\leftarrow$ PlanarSplitInterpolate(_Image_);

9 $\textnormal{{RefPlane}}\leftarrow\textnormal{{Planes}}[G]$;

10 $\textnormal{{WorkingPlanes}}\leftarrow\textnormal{{Planes}}[R,B]$;

11 foreach  _$\textnormal{{plane}}\in\textnormal{{Planes}}$_ do

12 EqualiseHistogram(_plane_);

13

14 end foreach

15foreach  _$\textnormal{{plane}}\in\textnormal{{WorkingPlanes}}$_ do

16 $\textnormal{{coeffs}}\leftarrow$ ($0.99$, $0.01$, $0.0$, $0.0$);

17 $\textnormal{{bounds}}\leftarrow$ ($1.0\pm 0.05$, $0.0\pm 0.05$, $0.0\pm
0.05$, $0.0\pm 0.05$);

18

19 repeat

20 $\textnormal{{distorted}}\leftarrow$ RadialDistort(_plane , coeffs,
Intrinsics_);

21 $\textnormal{{mask}}\leftarrow$ RadialDistort(_Blank(__) , coeffs,
Intrinsics_);

22 $\textnormal{{diff}}\leftarrow$ AbsoluteDifference(_distorted , RefPlane_);

23 $\textnormal{{error}}\leftarrow$ MaskedAverage(_diff , mask_);

24 $\textnormal{{coeffs}}\leftarrow$ BoundedBFGS(_coeffs , error, bounds_);

25

26 until  _StoppingCriteria(_error_)_;

27 end foreach

28WriteDatabase(_coeffs_);

29

30if  _Immediate_ then

31 WriteImage(_distorted , Metadata_);

32

33 end if

Algorithm 3.1 Correction coefficients recovery algorithm

The coefficient recovery algorithm estimates the correction coefficients which
give the minimum LCA, thus either artificial, e.g. calibration target or real
world images are suitable input data for this algorithm, given in Algorithm
3.1. This is broken down into discrete functional steps, shown in Figure 3.3.

start acquire source image decode image metadata load intrinsic lens
parameters histogram equalisation distort working plane compute mask compute
absolute difference compute average pixel value refine distortion coefficients
evaluate stopping criteria save data finishdescentconverged Figure 3.3:
Coefficient recovery steps

Acquire source image: An image acquired in raw format is selected and the data
loaded into a buffer. Lossy-compressed images (e.g. JFIF’s JPEG [38]) are
typically not used, due to two main reasons: primarily, non-raw formats
introduce colour aliasing during CFA demosaicing, due to adaptive algorithms
which adjust the local pixel colour due to directional gradients, or
misguidance artifacts, such as in AHD [42] — this prevents correct micro-
adjustment with the error minimisation step later. Secondly, compression
artifacts and chrominance subsampling will suppress high-frequency spectral
detail and attenuate low-amplitude edges, needed for optimal convergence. All
digital SLRs and most bridge and some high-end point and shoot cameras can
capture images in raw format, thus this is not restrictive for the scope of
this work (line 1).

Simple demosaicing is performed during this step to extract pixel value
information for each group of cells, performed using DCRAW111DCRAW is
ubiquitous for raw conversion integrated in a wide variety of open and closed-
source projects [21] to three separate R, G and B plane buffers as brightness
values, thus not scaled as luminance values according to the photopic equation
[18] (line 4).

Lastly in this step, linear interpolation is performed. Since the sampling
points at the three primary wavelengths are not co-incident but adjacent,
linear interpolation must be used to give continuous colour planes with the
same spatial size and position. CFA sensors have twice as many green-filtered
photosensitive cells than red or blue, thus the interpolated green plane can
have twice the spatial frequency than the red or blue planes. Modern
algorithms used to avoid colour ‘zipper’ effects at colour boundaries, such as
Hirakawa’s popular Adaptive Homogeneity-Directed demosaicing algorithm [42]
must be avoided, as these use colour selectivity to give homogeneous edge
features by detecting and minimising differences across colour planes (line
4).

Decode image metadata: Particular parameters specific to the lens settings at
the time the image are captured and encoded in EXIF tags. These are extracted,
and include the lens descriptor, focal length, aperture value and image
sequence number and are needed for associating the correction coefficients
with lens parameters for later use (line 2).

Load intrinsic lens parameters: Optionally if decentering has been computed or
is known, the camera lens decentering values, as part of the precomputed lens
intrinsic parameters are loaded from the database using the lens descriptor
(line 3).

Histogram equalisation: Camera systems seek to use the correct exposure to
achieve a wide global tonal range at image capture time; this is done in a way
to preserve the relationship between channels. Since the algorithm uses inter-
channel contrast, per-channel histogram equalisation is performed to ensure
convergence is stable on pictures with heavy colour casts, such as shot with a
non-white primary light source, e.g. under a sodium street lamp (line 8).

Distort working plane: Distortion of the working plane to the distortion
buffer is performed, using the current set of distortion coefficients.
Bilinear interpolation is used to minimise aliasing from non-integer pixel
displacements (e.g. that nearest neighbour would cause), and is selected as a
tradeoff of compute time and aliasing error (line 14).

Compute mask: A mask is computed to select only the subset of pixels common to
the reference plane and distorted plane; this is achieved by testing each
pixel in the distorted plane and setting the pixel in the mask if it is
interior to the distortion, otherwise clearing it. This later ensures that
pixels outside the resulting distorted frame do not influence the computed
pixel value average, ensuring the convergence will be the true minimum error
(line 15).

Compute absolute difference: For each pixel, the absolute difference between
the reference plane and distorted plane is calculated and stored (line 16).

Compute average pixel value: The non-masked absolute difference pixels are
summed and divided by the total non-masked pixel count in the distortion
buffer, giving the average per-pixel difference. This serves as a size-
independent measure of how homogeneous the planes are, required for
convergence (line 17).

Evaluate stopping criteria: The stopping criteria is evaluated to determine if
the rate of change of error is sufficiently low from the minimum observed, or
a reasonable limit on iteration count is reached. This is implemented as part
of the error minimisation algorithm selected, in the next step (line 19).

Refine distortion coefficients: The distortion coefficients are modified based
on feedback from the error function, via the L-BFGS-B method [13]. L-BFGS-B
seeks to minimise the value of the error function by adjusting the distortion
coefficients which are passed to the error function, which performs the per-
plane distortion and averaging of the inter-plane absolute difference to
compute the resulting error. L-BFGS-B was selected for faster and more robust
convergence than the Modified Powell method [9] and Nelder and Mead’s downhill
simplex method [67], due to the function being minimised having a low rate of
change over small changes in area [65]. Other gradient descent algorithms are
unsuitable due to inability to compute the function gradient or Hessian
matrix, which is inherently not possible from using an abstract error function
(line 18).

Save data: After convergence, selected image metadata, extracted from EXIF
fields in the image headers is stored along with the converged correction
coefficients for later use (line 21).

The coefficients given by minimisation are used to distort the original image
and written to a file, thus without the intermediate steps of histogram
equalisation. This part of the step is optional, and omitted when just
generating correction data for the database (line 23).

###  3.4 Image correction

The workflow of leveraging the correction database to correct pictures is
given in Figure 3.4 and broken down into discrete steps:

Decode lens parameters: Parameters specific to the image are decoded from the
image’s EXIF tags. This includes the lens descriptor, focal length and
aperture value. The database records matching the lens descriptor are
selected.

Search for nearest neighbour: For each record, a distance vector is computed
based on the difference between focal length, focus distance and discrete
aperture setting. This is part of the Nearest Neighbour Search [19] in the
next step, and can use cost vectors to optimise selection.

After iterating through all database records, the nearest neighbour is known,
and the correction coefficients are selected from this record.

Decode image content to RGB planes: The luminance values at each CFA cell are
separated into image planes by the colour filter pattern. Simple bilinear
interpolation is used to generate continuous colour planes with the same
spatial size and position. This step is achieved using DCRAW [21].

Distort plane: For each of the red and blue planes, the plane is resampled
into a new planar buffer using the selected correction coefficients. Accurate
resampling is achieved via a windowed approximation of the classic sinc
function [84], such as the Lanczos function [59]; a higher node count in the
windowed section can be used for more accuracy. Since the green plane acts as
a reference for other planes, and has a significantly higher spatial frequency
limit than the other planes, it does not undergo any transformation, thus
minimising loss due to aliasing and rounding to the colour depth.

Encode image data and tags: The planar image data is assembled into an RGB
triplet buffer and encoded into the preferred intermediate and lossless image
format. A suffix is added to the filename to associate the correct image with
the original and the image data is written. Additionally, the EXIF tags are
written and updated with information that this image has undergone LCA
correction applied, thus should not be reapplied.

start decode lens parameters search for nearest neighbour for all entries
decode image content to RGB planes distort plane for R,B planes encode image
data and tags finish Figure 3.4: Image correction with model database

###  3.5 Novel approach to correction quantification

Where an image is presented with a known (computed, predefined or otherwise)
distortion, the most direct measure of success after application of the
correction algorithm, is some quantified difference or _distance vector_
between the known distortion and the converged correction. A sensible metric
for this is the area of difference of the two polynomials representing the
actual distortion and the correction from the selected distortion
coefficients. In real-world situations, a more useful metric is the
_perception_ of chromatic aberration, since the ground-truth correction
coefficients are not known. Though perception cannot be accurately quantified,
the changes to image spatial frequency can be accurately quantified. Since the
HVS is sensitive to high spatial frequency in image luminosity, image clarity
or acuity is strongly linked to changes in spatial frequency [8]. If chromatic
distortion is introduced, per-plane spatial frequency will remain the same,
though the spatial frequency of the image’s luminosity will increase,
introducing false detail and reducing desired detail contrast.

Methods similar to this have not been encountered in this field, so a method
to quantify this change has been developed, optimised to this application.
Firstly, the image is transformed from a triple-plane RGB image into single-
plane luminance image via standard weightings from the CIE 1924 day photopic
equation, from Equation 2.9 [18]. The luminance plane is transformed into the
frequency domain through applying a Discrete Fourier Transform (DFT), giving
an array of complex numbers in the frequency domain:

$\displaystyle F_{1}(u,v)$ $\displaystyle=$
$\displaystyle\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(\frac{ux}{M}+\frac{vy}{N})}$
(3.2)

where $f(x,y)$ is the source image in $x$ and $y$ spatial domain with size $M$
and $N$ elements respectively. $F_{1}(u,y)$ is the destination image in $u$
and $v$ frequency domain. A second DFT is produced from the image being
compared against.

In order to generate a useful spectrogram, the algorithm proceeds to iterate
through the $u$ and $v$ dimensions in frequency space; it computes the spatial
frequency as the Euclidean distance from the 0Hz point and adds the magnitude
of the complex value at that point to the bin covering that spatial frequency,
producing a spectrogram $S$. Represented as a function of spatial frequency
$f$:

$\displaystyle S(f)$ $\displaystyle=$ $\displaystyle\sum_{\theta=0}^{\pi}\mid
F_{2}(f\cos\theta,f\sin\theta)\mid-\mid F_{1}(f\cos\theta,f\sin\theta)\mid$
(3.3)

where $S(f)$ is the spectrogram of spatial frequency $f$, and $F$ is the DFT
image.

Lastly, a windowed Hann filter is employed over the resulting one-dimensional
vector to attenuate undesired high-frequency components for the sake of
visualisation [40] and the result is plotted against radius, where the centre
is 0Hz, i.e. the whole-image average and up to the Nyquist frequency, thus one
pixel spatial resolution [20]. The resulting graph plots luminance spatial
frequency from the given input image. Where the input image is the subtraction
of two DFTs, this allows inspection of frequency changes introduced by LCA
occurring in the lens system, along with magnitude information. This method
for quantifying the loss in image acuity due to LCA is represented in
Algorithm 3.2.

1

2

3

input : Images to be analysed

output : Graph showing change in spatial frequency

4

5Ia $\leftarrow$ LoadImage ;

6 Ib $\leftarrow$ LoadImage ;

7 La $\leftarrow$ Luminosity(_Ia_);

8 Lb $\leftarrow$ Luminosity(_Ib_);

9 foreach  _image $\in\textnormal{{La}},\textnormal{{Lb}}$_ do

10 image $\leftarrow$ Pad(_image max (image.width, image.height)_);

11 complexMap $\leftarrow$ DFT(_image_);

12 foreach  _$x,y\in$ complexMap_ do

13 elem $\leftarrow$ complexMap $[x,y]$;

14 magnitude $\leftarrow$
sqrt(_${\textnormal{{elem}}_{i}}^{2}+{\textnormal{{elem}}_{j}}^{2}$_);

15 $\textnormal{{outputMap}}[\textnormal{{image}},x,y]\leftarrow$
log(_$\frac{1+\textnormal{{magnitude}}}{\textnormal{{image}}.\textnormal{{height}}}$_);

16

17 end foreach

18

19 centreDFT(_outputMap $[\textnormal{{image}}]$_);

20

21 end foreach

22diff $\leftarrow\textnormal{{outputMap}}[A]-\textnormal{{outputMap}}[B]$;

23 histogramMap $\leftarrow$ polarToCartesian(_diff_);

24 foreach  _$x\in$ histogramMap_ do

25 $\textnormal{{histogram}}[x]\leftarrow$ sumColumn(_histogramMap , $x$_);

26

27 end foreach

28filtered $\leftarrow$ filterHanning(_histogram_);

29 graph(_filtered_);

Algorithm 3.2 Quantification algorithm

###  3.6 Minimiser selection

The function measuring the difference between image planes for a given set of
distortion coefficients, shown in Equation 3.1 is minimised by varying the
coefficients used for distorting the image and measuring the error. Since the
function operates by taking the absolute difference of the image planes, it is
not possible to compute the function derivative or Hessian matrix, limiting
the choice of minimisation algorithms. Understanding the surface geometry of
this function gives insight into the complexity of the problem space and if
convergence is expected to be robust. To show the surface, brute-force
sweeping is conducted in the first and second order coefficients for the Red-
Green and Blue-Green plane pairs, shown in Figure 3.5. It can be seen that
there is a clear global minimum, thus robustness is expected within the
limits.

(a) red to green displacement (b) blue to green displacement

Figure 3.5: Error varying the first and second correction coefficients

A number of optimisation methods were evaluated for convergence performance
and stability. The simplex downhill algorithm [67] was found to not stabilise
on the global minima, but diverge on particular images. This was found to be
in part from the gradient descent also being designed for minimising linear
problems, whereas Figure 3.5 shows the non-linearity. Powell’s [77]
alternative to the simplex algorithm showed fewer iterations were needed for
convergence, however it exhibited stability issues with the noisier blue
channel with three pictures from the test set. The Broyden-Fletcher-Goldfarb-
Shanno (BFGS) algorithm [83, 34, 31, 12] was found to address these stability
issues due to taking smaller steps, limiting overshoot and non-linearity
error. A final refinement was to use the Limited-memory BFGS-Bounded
(L-BFGS-B) constrained optimiser [13].

Whilst careful selection of constraint values is important for consistent
stability, it was found that limiting the polynomial coefficients to a level
that would give in excess of maximum possible error (i.e. more than 30 pixels)
still maintained the desired stability; thus $1\pm 0.05$ was chosen for the
zeroth order coefficient and $|0.05|$ was chosen for higher order
coefficients; even the most basic of lens designs would have polynomial
coefficients within these bounds.

###  3.7 Summary

The methodology was broken into two phases, coefficient recovery and image
correction. Offline image analysis is performed in the coefficient recovery
stage; lens setting parameters are extracted and correction coefficients are
found and stored in the database. Image correction takes sets of image
parameters and selects the nearest match for the correction of unseen images
without incurring the coefficient recovery overheads, and runs the correction
algorithm.

Alternative approaches may involve discrete feature points, yielding far less
information and thus correction accuracy or stability, or may examine gray
areas. This approach performs no segmentation, using all areas, since non-gray
areas typically have colour components in two of the three red, green and blue
planes.

##  Chapter 4 Methodological evaluation

###  4.1 Validation

This chapter introduces a number of steps taken to validate the proposed
algorithm at levels of increasing complexity.

####  4.1.1 Chequerboard test image

Initially, a classic chequerboard image is rendered; this is distorted via
radial pixel remapping from Equation LABEL:eqn:seidel-approx. Coefficients
presenting a significant level of LCA distortion are selected:

$\displaystyle r_{dest}$ $\displaystyle=$ $\displaystyle
a\times{r_{src}}+b\times{r_{src}}^{2}+c\times{r_{src}}^{3}+d\times{r_{src}}^{4}$
$\displaystyle a$ $\displaystyle=$ $\displaystyle 0.98$ $\displaystyle b$
$\displaystyle=$ $\displaystyle 0.01$ $\displaystyle c$ $\displaystyle=$
$\displaystyle-0.01$ $\displaystyle d$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\therefore\ \ r_{dest}$ $\displaystyle=$ $\displaystyle
0.98r_{src}+0.01{r_{src}}^{2}-0.01{r_{src}}^{3}$ (4.1)

The correction algorithm is executed on the resulting distorted image; it is
expected that it will converge on the minima at which there is least LCA
distortion, i.e. back to LCA-free. This can only be tested with no uncertainty
because the pattern is achromatic and by definition, LCA-free to start with.

Crops of the test image are shown in Figure 4.1. By inspection of the
righthand column of crops, no residual LCA artifacts, which would appear as
linear colour fringes, are visible down to sub-pixel level. This shows three
results; firstly the algorithm has converged on the global minima and not at
any local minima. Secondly, convergence is accurate due to correction at the
sub-pixel. Lastly, convergence is reliable across various radii as no LCA is
detectable in any of the resulting crops.

Figure 4.1: Cropped sections from original (left) distorted (center) and
corrected (right) test image

Visual inspection is not reliable, due to a number of factors: firstly,
subjectivity and perception may label measured error insignificant. Secondly
for the error to be observed, it needs to be perceptually and psychovisually
significant, i.e. colour offset within a noisy texture generates a low HVS
response compared to high-contrast edges, thus is harder to visually flag.
Lastly, quantifying the magnitude of the change is important to understand and
measure the performance of the algorithm. To address these concerns, Algorithm
3.2 is employed, measuring the change in spatial frequency from the original
aberration-free pattern relative to the distorted pattern. The quantification
algorithm produces a graph with the magnitude of the change in spatial
frequency plotted against spatial frequency.

The process of artificially introducing effectively-LCA distortion, seen as
the green line in Figure 4.2, increases low spatial frequency up to around
0.15 Nyquist, and decreases the higher frequency magnitude significantly up to
around 0.36 Nyquist. Since the chequerboard image has regular intervals and
maximal contrast at transitions, there is a very strong response at harmonic
frequencies of the block size, due to the periodic nature of DFTs. In
essentially damping the luminance contrast and geometrically spreading the
contrast transitions, attenuation measured in the DFT at higher frequencies
would be expected and is observed. The distortion recovery algorithm, shown as
the blue line, effectively reverses the spatial frequency changes, including
restoring the higher frequencies that were attenuated in order to converge on
minimal distortion. From the graph, at $f=0.18N$ where $N$ is the Nyquist
frequency of the image, we see an around 20% reduction in frequency
restoration. This is attributed to aliasing loss, since the pattern observed
is symmetrical with respect to the introduced distortion.

Figure 4.2: Change in spatial frequency from undistorted to distorted (blue)
and distorted to corrected (green) pattern images

####  4.1.2 Ray-traced scene

As a second validation step, a ray-traced scene representative of real-world
luminance, spatial detail and colour complexity is selected and rendered,
shown in Figure 4.3. This provides a chromatic aberration-free source image,
which is consequently distorted in the same manner with known coefficients to
test the performance of the algorithm with typical spatial frequencies
present. Since the image is not acquired from a Bayer array, it has higher
information density and approaches the Nyquist limit. As a result, more detail
loss due to interpolator quality, aliasing, and approximation or accuracy
limits in the image warping is expected.

Figure 4.3: Rendered test image

Figure 4.4: Cropped sections from original (left) distorted (center) and
corrected (right) rendering

Distortion on the ray-traced scene is executed with the known coefficients,
then the correction algorithm is iterated, shown in Figure 4.4. Figure 4.5
illustrates the increase in luminance spatial frequency due to LCA, shown in
blue. The high frequencies approaching half-Nyquist (i.e. one-pixel width)
have been attenuated as expected. The algorithm’s correction, shown in green
illustrates excellent performance at lower spatial frequencies, mirroring the
shape and magnitude of the distortion with approaching zero frequency loss. At
spatial frequencies above $f=0.2N$, significant attenuation to the correction
is seen. Since all the pixels in the image have been remapped with
interpolation and antialias filtering, this is expected. Tests conducted in
the next section will identify the impact of this on real-world image data,
where spatial frequency from a lens system and Bayer array is far less than
half-Nyquist.

Figure 4.5: Change in spatial frequency from undistorted to distorted (blue)
and distorted to corrected (green) rendered images

###  4.2 Performance analysis

####  4.2.1 Real-world test image

Firstly, the real-world example previously shown in Figure 1.4 to illustrate
LCA was selected. The LCA-correction algorithm was applied, shown in Figure
4.6. No residual LCA can be visually observed from the crops, which is the
desired outcome. Since simple linear interpolation was used and not more
advanced demosaicing algorithms, expected ‘zipper’-like artifacts are
observed; for analysis without this, pixel remapping would have to be
performed prior to demosaicing. If this were performed, advanced demosaicing
would enhance quality further, rather than simple linear demosaicing.

Figure 4.6: Cropped sections from original (left) and corrected (right)
photograph

To visualise the changes, difference maps are generated by taking the absolute
difference of each corresponding pixel in the blue and green planes, and
separately for the red and green planes; white edges would be seen where
chromatic differences exist including due to LCA. It is worthwhile noting that
even for an image with zero LCA (e.g. an artificial image), there will be a
constant difference map total area which the correction algorithm would not be
able to minimise further, thus there is a fixed offset in addition to the
edges due to the LCA. As a result, it is seen from Figure 4.10 that there is a
significant reduction in structure comparing the uncorrected to corrected
difference maps, particularly in the blue-green plane difference.

Next the quantification algorithm was executed to quantify the spatial
frequency changes between the original and corrected image, shown in Figure
4.7. A significant reduction in spatial frequency of peak 0.015dB is observed.
This translates into information which is more useful to the HVS both in
luminosity and chromaticity, and the increase in homogeneity can be observed
by the picture looking ‘clearer’ as a result.

Figure 4.7: Change in spatial frequency due to LCA correction

####  4.2.2 Run-time cost

Since the correction and quantification algorithms are highly computationally
intensive and intentionally avoid subsetting the data they measure for
robustness, the runtime can be consequently significant, shown in Figure 4.8;
time was measured on an Intel Core i5-661 desktop platform.

For the correction algorithm, the runtime cost is dominated by the image size
(in pixels) and the efficacy of the L-BFGS-B algorithm minimising the error in
a non-linear multivariate problem; this influences the number of iterations of
the relatively expensive error function over the image. The majority of time
within the error function is spent on non-linear pixel remapping, which
produces a lot of out-of-cache access. OpenCV is internally used for the
remapping without a highly-optimising compiler, thus this could benefit
significantly from optimal load/store scheduling and ideal prefetch scheduling
to avoid frequent processor pipeline stalls due to waiting for data to arrive
from the memory controller. The left-most and small blue bar indicates image
loading and plane splitting (0.7%); the green bar shows the total time
subtracting image planes and computing the masked average pixel value (5.4%);
the red bar represents the time generating the pixel remapping offsets
performing the interpolated lookups (93.9%).

The quantification algorithm has a far shorter runtime, as no iteration is
performed; the left-most and small blue bar shows that loading both input
images has a similar overhead as in the correction algorithm (5.3%). The two
DFT transformations shown in the green bar account for 15% of the overall
time, while the rest of the computation time is shown in the red bar (79.7%)
and is due to binning and accumulating the DFT output based on radius, i.e.
spatial frequency.

5010015020025025737.80secondscorrectionquantification Figure 4.8: Runtime
breakdown of algorithms

####  4.2.3 Suite of images

The improvement is reflected in the reduction in the area of the plane-
difference maps, compared in graphical form in Figure 4.9. The reduction in
the red to green and the blue to green difference map area are as shown in
Table LABEL:tab:edge-reduction. This illustrates just the reduction in
difference areas the correction algorithm was able to achieve; inspection of
the resulting corrected images confirmed that all the minima were true global
minima and not leading to misconvergence, which would manifest as introduction
of additional LCA artifacts. Since the difference area without LCA is not
known, this relative reduction in difference area carries no absolute meaning.
Later, the quantification algorithm will be employed to measure differences
with more significance to the HVS.

0510152025303540457.98% ave3.86% ave12.1%
ave1234567891011121314151617181920212223242526test image indexpercentage
reduction Figure 4.9: Difference map area reduction. Blue bars show blue-green
edge reduction; red bars show red-green edge reduction.

The correction algorithm was executed for all of the images in the bank of
representative real-world test images, shown earlier in Figure 2.6. The change
in spatial frequency for all images was quantified and plotted together, shown
in Figure 4.11.

[botcap,caption=Difference map area reduction,label=tab:edge-
reduction]lllImage & R-G reduction B-G reduction 1 7.33% 12.5% 2 6.60% 12.3% 3
4.45% 6.95% 4 7.14% 19.0% 5 6.49% 21.7% 6 3.63% 10.3% 7 10.9% 12.6% 8 13.7%
29.2% 9 6.00% 9.05% 10 2.30% 10.8% 11 3.26% 12.3% 12 7.73% 11.0% 13 0.31%
5.59% 14 0.21% 2.32% 15 1.88% 7.79% 16 0.80% 5.63% 17 0.38% 8.84% 18 3.37%
9.17% 19 1.29% 26.4% 20 1.08% 11.0% 21 2.18% 12.4% 22 1.82% 14.6% 23 1.49%
7.85% 24 3.90% 15.4% 25 1.64% 12.0% 26 0.567% 7.56% Average 3.86% 12.1%

Figure 4.10: Cropped sections from uncorrected red-green (far left), corrected
red-green (left), uncorrected blue-green (right) and corrected blue-green (far
right) difference maps Figure 4.11: Change in spatial frequency for test
images. Colour shows image index.

###  4.3 Comparison with existing methods

Three images were selected for testing against a widely-used solution for LCA
elimination: Adobe Photoshop Creative Suite 4 [45]. The demosaiced images that
were used for proposed correction were also used in the Photoshop LCA
correction to avoid comparing differing raw converter and internal processing
steps. The procedure employed iteratively refined the three ‘Green/Magenta
Fringe’, ‘Blue/Yellow Fringe’ and ‘Red/Cyan Fringe’ controls, until LCA
artifacts were eliminated as far as possible by inspection of a 200%
magnification of the image at a range of places throughout the image.
Following this, quantification was employed to measure the improvement. The
results are analysed in the next section.

Figure 4.12: Cropped sections from original (left) Photoshop (center) and
proposed correction (right) Figure 4.13: Change in spatial frequency for
Photoshop (blue) and proposed (green) correction in the first selected image

Figure 4.14: Cropped sections from original (left) Photoshop (center) and
proposed correction (right) Figure 4.15: Change in spatial frequency for
Photoshop (blue) and proposed (green) correction in second selected image

Figure 4.16: Cropped sections from original (left) Photoshop (center) and
proposed correction (right) Figure 4.17: Change in spatial frequency for
Photoshop (blue) and proposed (green) correction in the third selected image

###  4.4 Results discussion

By inspection of Figure 4.11, there are a number of characteristics observed,
given by:

  * •

Spatial frequency is observed to tend towards zero at 0 Hz and half-Nyquist.
If this were observed to not be the case, error would have been introduced
from the correction algorithm, or mismeasurements from the quantification
algorithm

  * •

There is no spatial frequency (e.g. texture complexity) increase, i.e. all the
changes serve to reduce spatial frequency; this corroborates with the basis of
chromatic regions gaining homogeneity and minimising additional spatial
frequency due to the reduction of the resultant aliased artifacts

  * •

There are two observable peak reductions around $0.12N$ and $0.35N$, where $N$
is half the sensor Nyquist frequency; these would correspond to LCA features
around 1.4 pixels and 4.1 pixels in width, which fits with the LCA feature
size of around 4 pixels seen in Section 2

  * •

For this lens system, Figure 4.9 shows that significantly more blue edge area
is minimised than red edge area from Figure 4.9; this suggests that the lens
system is optimised for lower LCA at redish wavelengths around 700 nm than
blueish around 470 nm, which is a plausible design optimisation for the
general use-profile of this lens. An areal reduction on average of 7.45% was
achieved

Another technique was initially evaluated and proven to give poor robustness
for test image corrections; this algorithm generated a list of image features
with Lowe’s SIFT [61] for each image plane, discarding outliers and employed
feature correspondence to give an error measure, subsequently minimised via
Least-Squares [64] to converge the image planes. It was found that there was
an insufficient number of useful image features to attain complete and
accurate convergence in a significant number of the test cases; while others
did not converge to the global minima, but to a local minima, presenting
miscorrection. It was found that detecting discrete image features may not
give coverage at a sufficient number of different radii from the image centre,
which would optimise the image at a sub-range at expense of under-measured
areas. The proposed technique avoids these issues by minimising the continuous
planar difference, so each differing pixel contributes to the error measure.
This provides a far more stable and predictable descent via the error
function. Discrete keypoints were used by Cecchetto [14], however the use of
statistical correlation keypoint pruning was shown to introduce a level of
uncertainty leading to subsequent correction error.

From the bank of test images, the proposed algorithm exhibited robust
behaviour and showed no misconvergence. This is due to the distortion having
radial homogeneity, making the approach sensitive to minima from features
which are perpendicular to the line towards the centre of the image. As the
lateral chromatic aberration is an image-global phenomena, after taking into
account any lens decentering, it is experienced uniformly at all angles at a
given distance from the centre of the image, thus local LCA features
contribute small influence over the error minimisation surface. This property
is fully exploited by performing the algorithm over the whole image, rather
than just a subsector of the image, or on discrete image features. Whilst it
is entirely possible to prepare a theoretical image which would defeat this
mechanism e.g. via red, green and blue concentric circles of differing but
near radii, this has not occured in practice, or nor is likely to. Per-channel
histogram equalisation was introduced as a step in the algorithm to balance
out the potentially uneven influence each channel e.g. due to a colour cast
from a non-white light source.

For extremal or unexpected cases, the algorithm limits coefficient excursions
to plausible values by employing bounded error minimisation [106]. For any of
the distortion polynomial components to be greater in magnitude than 0.05,
creates a level of displacement that would only be useful in correcting an
unrealistic level of LCA, thus limits are retained. Secondarily, this provides
a set of constraints for the optimiser to work within, allowing searching for
the minima to converge quicker than an unconstrained search [106].

At present, the demosaicing step is prior to the displacement correction,
itself which uses linear interpolation to avoid introducing additional
uncertainty and artifacts. From the uncorrected and corrected crops zipper-
like artifacts are consequently visible at chromatic edges arising from the
simple linear interpolation not exploiting spatial locality to detect local
contours, as done in more advanced demosaicing algorithms, as previously
addressed. With the impact of LCA present in a picture, at chromatic edges,
advanced demosaicing would introduce false edges due to the demosaicer design
goal of local homogeneity; should the correction step be performed prior to
demosaicing, even further gains would be achieved.

###  4.5 Summary

Validation of the proposed LCA minimisation algorithm was performed
incrementally with an artificial checkerboard pattern, followed by an
artificial ray-traced image to prove basic operation. These were used to
provide the basis for real-world testing, and provide an initial model within
which the real-world examples would correlate. Following this, the initial
example correction was presented, followed by correction of the test images.

The data from the image correction was quantified and graphed, allowing
numerical inspection of the results. The data on the graphs show that higher-
order luminance spatial frequencies were attenuated, indicating false aliasing
from LCA was removed. Visual inspection confirmed this, and numeric comparison
with a leading method showed a significant improvement with this novel
approach.

##  Chapter 5 Conclusion and further work

###  5.1 Review

In this work, the imaging performance of practical photography systems was
reviewed and current limitations were examined. The poorly addressed and
increasingly significant limitations that LCA introduce were analysed in
detail. Various solutions have been developed [101, 99, 53, 23, 63, 93, 57,
62, 32, 78, 60, 14], however many present limitations on generality,
application area or impose impractical constraints for the correction of
arbitrary images.

Real-world examples of LCA in photographs were presented along with discussion
of parameter coverage. Next, the goals and scope of the solution needed were
defined. The correction algorithm methodology was then broken down and
presented into steps with justification.

A quantification framework was developed to measure LCA, then later the
correction algorithm was incrementally validated with test subjects from a
classic chequerboard test pattern and a high-complexity rendered scene. The
quantification framework was used to measure the improvement to artificial and
real-world test cases, and the performance of a leading LCA correction method
was compared with the proposed method. Results were presented that showed
visual elimination of LCA. Finally, the results including change in edge
difference and spatial frequency were presented and discussed.

####  5.1.1 Mechanics

The algorithm for the correction of LCA presented in this work is based around
minimising the sum of absolute difference between the reference green plane
and the target red and blue planes; this achieves the goal of minimising LCA.

Firstly, preprocessing the source planes with histogram equalisation was
performed to ensure the magnitude of the per-plane pixel values was
comparable; this was introduced to address the problem of converging on false
minima from images with strong colour casts, such as shot under sodium or
other poor artificial lighting.

Convergence on minimal LCA is then achieved by iteratively varying the Taylor
series coefficients in a radial distortion function; this is used to warp the
target plane to align up with the reference plane, mitigating the lateral
pixel displacement caused by differing paths due to diffraction. Error is
measured by taking the average absolute difference of the target and reference
planes for each pixel, with an inner boundary of the two planes to mask
outlying pixels. L-BFGS-B descent was employed to take advantage of a priori
knowledge of extreme LCA limits that would not be seen in practise. This
presented robust convergence and required fewer descent steps than other
techniques, which were less suitable due to the non-linearity of this problem.

Figure 5.1: SIFT used to detect feature points per channel

Notably, early testing showed that using discrete feature points failed to
provide adequate information for stable and accurate convergence on minimal
LCA; Figure 5.1 is an upper-right crop showing good feature-point correlation
in the shadowed area under the roof, but disparity elsewhere. The primary
advantage of the final algorithm is therefore the maximal use of image
information to achieve accuracy and robustness through calculating inter-plane
absolute difference. This also mitigated uncertainty and instability in the
decent process.

####  5.1.2 Results

The results are separated into three areas. Firstly, from setting the context
and presenting the algorithm steps and rationale, initial absolute difference
crops given in Figure 4.10 show the design principle — comparing the absolute
difference in two planes on the LCA-uncorrected image and the LCA-corrected
image. It can be seen that the corrected absolute difference map is overall
darker, due to more coincident subtraction; this concurs with expectation.

Next, a simple testcase was constructed by rendering a classic chess-board
calibration target. This was distorted with a Taylor-series polynomial and
corrected; crops were compared to validate the algorithm. Visually, the chess-
board converged to within sub-pixel expectations. Quantification in Figure 4.2
showed correct inversion of the frequency changes, though some frequency loss
was noted. Executing validation on a high-complexity ray-traced scene with
spatial frequency exceeding Bayer array capability showed that convergence was
visually accurate, such that no error could be observed in the crops. This was
understood due to the increased level of available visual information, and
highlights the architectural benefit of using maximum information from all
pixels. Quantifying the frequency change in Figure 4.5, there is notable
frequency loss due to aliasing error. Both observations of spatial frequency
loss are not significant when dealing with images from Bayer arrays; the red
and blue channels are sub-sampled by a factor of four within Bayer
architecture. The green channel is passed through unmodified, so experiences
no change.

Subsequently, LCA correction was performed on all of the images in the test
set given in Figure 2.6. Overall inspection confirmed no divergent correction.
As a result of misconvergence being observed in the twelfth test image due to
coloured glass windows while developing the algorithm, analysis revealed
descent to a false minima; the cause was found to be the uneven relative
weighting of the image planes being differenced. Addition of the per-channel
histogram equalisation step addressed this. A further step was taken to
introduce non-linear bounded minimisation via L-BFGS-B minimisation algorithm.

Inspecting the corrected test images at pixel-level confirmed that convergence
was accurate such that the visible chromatic error was wholly due to the
simple (linear interpolation) demosaicing algorithm; this demonstrated that
the minima was central within error from linear interpolation demosaicing.
Quantification of all the images in Figure 4.11 showed an interesting trend:
that there are two peaks in spatial frequency reduction, at around 0.15
Nyquist and 0.35 Nyquist. This is understood to be due to the frequency and
geometry of LCA feature size. Nevertheless, all pictures responded positively
to improvement.

In order to test algorithm performance against a well understood existing
approach, LCA correction of the images was undertaken in Photoshop CS4; best-
effort manual correction was executed, and the same source image was processed
with the proposed algorithm. Crops were extracted and compared. Visually, LCA
is shown to be eliminated as far as possible in both corrected images, as
expected. Employing quantification revealed a 20% to 22% peak reduction in
spatial frequency over the Photoshop approach, thus converging on chromatic
homogeneity down to the sub-pixel level in an optimal manner.

####  5.1.3 Summary

A number of significant observations and findings were drawn from the outcome
of this work:

  * •

When correcting a large sequence of pictures, it was found that there were
exactly three lens parameters at capture time that affected the LCA correction
coefficients: focal plane, focus distance and aperture size. The first
parameter is natural, as there is movement of lens elements when this is
varied. The last two parameters are less intuitive when varied, since focusing
moves internal lens elements a short distance; changing the aperture size
doesn’t move lens elements. Critically, both do alter the path of light.

  * •

Detecting discrete feature points on the image colour planes and finding
correspondence from the feature descriptor similarity was found to locate the
order of hundreds of inter-plane stable features, and was used in the initial
approach. Later this proved to be inadequate for accurate convergence on the
minimum-LCA solution. Following this, the architecture was modified to exploit
areal differences to deliver the sensitivity and accuracy anticipated.

  * •

Performing LCA minimisation on a subsampled image (e.g. quarter-size)
demonstrated the same robust convergence on the true minimum-LCA solution,
though an increase in error at a sub-pixel level was observed. Since
performance was secondary to accuracy and correctness in this work, this was
left aside as a future potential optimisation, using phased convergence with
various pyramid sub-sample levels.

  * •

A number of other camera systems were tested during late development with only
a small set of images; as expected, LCA correction showed correct convergence
and produced satisfying results

  * •

The algorithm proved to deliver robust convergence on minimum LCA, while still
allowing very wide convergence bounds for experimental lens systems

  * •

Convergence was accurate to the sub-pixel level to the extent that demosaicing
was the limiting factor for image quality; this limitation could be surpassed
by performing LCA minimisation and correction prior to demosaicing

  * •

From the set of test images shown in Figure 2.6, 3.86% average edge reduction
was found for the red plane and 12.1% for the blue plane; this indicated the
Nikkor lens used to capture the test images is optimised for lower red-green
LCA than blue-green LCA

  * •

Comparable attenuation was observed at higher spatial frequencies in the
proposed method; further analysis with a rendered target and Bayer simulation
would allow accurate quantification of this loss, and a similar level was
observed from the Photoshop method

  * •

Finally, significant image noise was found to not disrupt or hinder accurate
convergence; this is considered due to the stochastic nature of sensor noise,
averaging to zero across small patch areas, small with respect to the whole
image area used in minimisation

####  5.1.4 Limitations

Due to the iterative and whole image analysis and warping, the cost of each
step is high, and thus the overall minimisation of LCA on an image is
relatively expensive: with single-threaded execution, it took between 100 and
300 seconds to process an image. For this reason, the framework to re-use
correction data was developed. It was noted that around 80% of processing time
(Pareto principle) is spent in the spatial image remapping; since there is no
data dependency, this is highly suitable for processing on a highly-parallel
system e.g. a Graphics Processing Unit (GPU), and would achieve almost linear
speedup.

This computational constraint precludes application onto an embedded hardware
target, such as a camera, unless implemented in a custom image processing
engine in silicon or with a Field Programmable Gate Array (FPGA). This doesn’t
present a problem in practise, since the lenses available for a particular
camera platform can be profiled and a database built, or correction performed
offline.

Lastly, due to spatially remapping the red and blue planes, error is
introduced from pixel aliasing; this was previously demonstrated with an
artificial rendered image. While Bayer arrays are used in the vast majority of
camera systems, some exist which capture full RGB data per-pixel, such as
using the Foveon X3 sensor. If a lens was employed which had a Nyquist
frequency exceeding the sensor’s Nyquist, correction would incur spatial
frequency loss in the red and blue planes for areas in focus. This is
considered a reasonable tradeoff, and is the cost of post-processing (i.e.
non-optical) correction. This limitation was also seen in the Photoshop
method.

###  5.2 Claims and contributions

The contributions of this work are divided into three areas:

  * •

An algorithm delivering optimal LCA correction via robust convergence using
arbitrary images. The associated error is minimised by the design of using all
the information available from the image at every iteration, ensuring accuracy
is maximised. The initial application of this work is photogrammetry, however
this can equally apply to the areas found in broader literature, notably
including microscopy and telescopy.

  * •

Mechanism to quantify the shift in spatial frequencies in a way useful to the
HVS was developed. This provided a platform to develop and test the LCA
correction algorithm, and allowed correction stability to be understood at low
frequency and high up to half-Nyquist spatial frequency, where visual
inspection would be clearly deficient.

  * •

A system for resulting correction data being associated with extracted image
parameters for efficient offline correction was presented. Realising this in
practise was enabled by locating near correction parameters, with further
improvement possible via multi-variate interpolation.

###  5.3 Future work

####  5.3.1 Enhanced interpolation

Though the presented algorithm uses Nearest Neighbour Search to locate the
most suitable of correction coefficients, this is clearly non-optimal for two
cases. Firstly, where the number of seen images and thus spread of correction
coefficients is low, the average distance between lens parameter values is
unfavourably high; this will correlate to higher average error, and may limit
whole-system correction accuracy. Error from lens element position hysteresis
and sensor discretisation are assumed small. Secondly, to achieve optimal
correction, some level of coefficient interpolation must be used. There are a
set of problems associated with this, which are briefly discussed: since the
data is not at regular grid intervals, it is considered sparse and thus
classic linear interpolation between data points cannot be used. Instead,
Delaunay triangulation must be performed in the three dimensions of the
parameter space to locate the tetrahedron covering the location in which the
point of interest lies [24]; this is the basis for Natural Neighbour
interpolation [85]. Further, polynomial interpolation must be used, as the
correction coefficients cannot be interpolated individually. Methods such as
Aitken-Neville [1, 68] would be suitable.

Lastly, a simple step to achieve less information loss from pixel aliasing, is
to use a more advanced interpolation algorithm for pixel remapping during
warping, such as in the sinc family, e.g. Lanczos [59].

####  5.3.2 Parameter modelling

It should be noted that as the error in the correction coefficients is low due
to accurate convergence, interpolation between points is sufficient rather
than construction of a model of the correction coefficients mapping onto the
lens parameters. Algorithm 5.1 is presented as a solution using first order
linear interpolation.

1

2

3

input : source image for LCA correction, complete with EXIF data

input : database of previously computed correction coefficients for a range of
images

output : interpolated correction coefficients

4

5Image $\leftarrow$ LoadImage(__) ;

6 LensIdentifier, Aperture, FocalLength, FocusDistance $\leftarrow$ Image.EXIF
;

7 Database $\leftarrow$ OpenDatabase(_LensIdentifier_);

8 Radius $\leftarrow$ Database(_“MaxRadius”_);

// graph average of both polynomials

9 for  _parameter
$\in\\{``Aperture^{\prime\prime},``FocalLength^{\prime\prime},``FocusDistance^{\prime\prime}\\}$_
do

// locate two nearest neighbours in Euclidean space

10 lowerDist $\leftarrow$ NearestNeighbour(_Database , parameter, lower_);

11 upperDist $\leftarrow$ NearestNeighbour(_Database , parameter, upper_);

// graph average of both polynomials

12 for  _interval $\leftarrow 0\ \textnormal{{to}}$ Radius_ do

13 a $\leftarrow$ lowerDist(_interval_);

14 b $\leftarrow$ upperDist(_interval_);

15 average [parameter ]
$\leftarrow\frac{\textnormal{{a}}\times(1-\textnormal{{lowerDist}})+\textnormal{{b}}\times(1-\textnormal{{upperDist}})}{2}$;

16

17 end for

18

19 end for

20

// use standard polynomial curve fitting

21 newCoeffs $\leftarrow$ PolyFit(_points $\in$ average_);

return(_newCoeffs_)

Algorithm 5.1 First-order coefficient interpolation

One solution for further consideration could be selected at the expense of
implementation complexity and may be more suitable if convergence was
inaccurate: parametric surface fitting to model the lens parameters. Earlier
steps create database entries for a particular lens, storing the computed
correction coefficients and relevant lens parameters at image-capture time.
The correction coefficients are computed by minimising an error function which
sums the vector of the fitted points to the ground truth data. Together, these
are used to derive the correction coefficients for an unseen image. The full
list of steps are:

Database storage: A database record is constructed with sufficient information
to identify the camera lens and related information such as lens mounting
(thus crop factor). This is later used to select the relevant correction
parameters for a particular lens and mount.

Parametric surface fitting: Once the data from each image in the input library
has been stored into the database, each database record is accessed to
identify what lens it was taken with; only images with the selected lens are
considered in the model, thus images taken with other lenses are not
considered and skipped.

A parametric surface equation is used to model the focal length and distortion
relationship. For a given focal length, the experimental approach Fang et al
[28] take, finds first-order fitting sufficient for their application; the
optical design theory calculated by Gross et al [37] concurs with this.
However it has been seen in practice, that simple first-order LCA correction
is insufficient. In order to design in scope for this, a higher-order
polynomial is used in both axes. Where second order is sufficient, the third
order coefficients will be observed to be small.

The resulting parametric model with the constant term and cross-products
removed gives six coefficients:

$z=bx+cy+dx^{2}+ey^{2}+dx^{3}+ey^{3}$ (5.1)

Marquardt’s least-squares fitting [64] is used to converge the coefficients,
minimising error.

Profile generation: The parametric surface coefficients are stored along with
the lens data, so can later be retrieved for offline correction.

####  5.3.3 Implementation performance

The last identified area for future work, is that an implementation leveraging
GPU hardware would enable interactive. As previously revealed, around 80% of
the processing time is spent in the plane remapping function, of which is
highly suited to GPU and highly-parallel architectures. Since there is no data
dependency between the microkernels that process patches of destination pixels
in the remapping function, the parallel scheduling and execution will
efficiently use processing resources. With fine-grained hardware scheduling,
throughput and efficiency far higher than on a general-purpose processor would
be seen, resulting in low execution time.

###  5.4 Summary and concluding remarks

LCA correction is clearly a powerful and compelling step to reduce the loss
that lens systems inevitably impose. It has been found that the proposed
solution addresses the need of general high-fidelity photography and provides
a solid foundation for further development, e.g. to realise faster processing
or correction parameter interpolation. Since the presented algorithm can be
applied in general to any imaging system, it has wide application across a
variety of fields, including but not limited to the significant areas of
telescope and microscope imaging.

The author believes in the free availability and openness of software and will
make the source code available for use, along with this thesis. A future step
would be to integrate this work with Zabolotny’s Lensfun optical correction
library [103], which itself is used in a number of popular image manipulation
tools for converting from raw sensor data. Photographers and scientists would
be able to generate and further submit high-quality lens LCA correction data
for the library to automatically use, further realising the benefit.
$\blacksquare$

## References

  * Aitken [1932] A. C. Aitken.  On interpolation by iteration of proportional parts without the use of differences.  _Proceedings of Edinburgh Mathematical Society_ , 3:56–76, 1932.
  * Alda [2003] J. Alda.  _Encylopedia of Optical Engineering_.  CRC Press, second edition, 2003.
  * amd Joaquín Otón and Pérez-Cabré [2006] M. S. M. amd Joaquín Otón and E. Pérez-Cabré.  Dynamic compensation of chromatic aberration in a programmable diffractive lens.  _Optics Express_ , 14(20):9103–9112, 2006.
  * Anagnostopoulos and Lesser [1995] C. N. Anagnostopoulos and M. P. Lesser.  _Cameras and systems for electronic photography and scientific imaging_.  SPIE, 1995.
  * Bayer [1976] B. E. Bayer.  Color imaging array. US patent 3971 065, 1976.
  * Benhimane and Mali [2004] S. Benhimane and E. Mali.  Self-calibration of the distortion of a zooming camera by matching points at different resolutions.  _Intelligent Robots and Systems_ , 2004.
  * Bockaert [2006] V. Bockaert.  Determination of chromatic aberration correction parameters, 2006.  URL http://www.dpreview.com/reviews/AdobePhotoshopCS2/page2.asp.
  * Bovik [2005] A. C. Bovik.  _Handbook of Image and Video Processing_.  Academic Press, 2005.
  * Brent [1973] R. P. Brent.  _Algorithms for Minimization without Derivatives_.  Prentice-Hall, 1973.
  * Brown [1966] D. C. Brown.  Decentering distortion of lenses.  _Photogrammetric Engineering_ , 32(3):444–462, 1966.
  * Brown and Lowe [2003] M. Brown and D. G. Lowe.  Recognising panoramas.  _International Conference of Computer Vision_ , 2:1218–1227, 2003.
  * Broyden [1970] C. G. Broyden.  The convergence of a class of double-rank minimization algorithms.  _IMA Journal of Applied Mathematics_ , 6(1):76–90, 1970.