Quaternionic phase and axis/colour retrieval

We demonstrate a quaternion-based phase and axis retrieval algorithm and a quaternionic realization of the Gerchberg–Saxton method, particularly suitable for RGB images. The RGB image is represented here as quaternion-valued image in polar form with the components quaternionic magnitude, axis and phase. We discuss their importance and relations in spatial and Fourier domain. We illustrate three versions of retrieval: the quaternionic phase only, the axis only and the joint phase and axis retrieval algorithm derived from the classical Gerchberg–Saxton algorithm. We discuss how the retrieval results are influenced by the choice of the Fourier axis of the quaternionic Fourier transform. Finally, we show that by including more constraints and prior knowledge, respectively, the retrieval can be improved.


Introduction
The Gerchberg-Saxton (GS) algorithm [1] is a well known algorithm to reconstruct the phase of an object function while only the magnitudes in spatial and Fourier domain are given. Such phase retrieval problems arise in many fields and found modifications over time [2]. In the case of complex grayscale images the GS algorithm typically converges very fast towards a good reconstruction of the image. Phase retrieval is also of interest for colour images such as RGB images. However, the algorithm has to be executed three times, i.e. for each colour channel of the image. While quaternion-based algorithms and the theory of quaternionic phase retrieval have been considered before [3,4], we present a new approach for quaternionic phase retrieval based on the B Martha L. Zimmermann martha-lina.zimmermann@math.tu-freiberg.de Swanhild Bernstein Swanhild.Bernstein@math.tu-freiberg.de Bettina Heise Bettina.Heise@jku.at 1 GS algorithm. This approach uses the quaternion polar form (2) and avoids splitting the RGB image into the three colour channels. The quaternionic representation allows the consideration of the colour image at once while only knowing the quaternionic magnitudes of the entire colour image in spatial and Fourier domain (FD). The use of quaternionic Fourier transforms (QFT) (8), (9) allows a modification of the classical (complex-valued) GS algorithm to a quaternionic version and herewith to be extended to colour image retrieval.
In the second section we will recall some of the fundamentals of quaternions and the QFT. The third section explains the representation of RGB values as quaternions. Furthermore, we will discuss the importance of magnitude, axis and phase of quaternions and quaternionic Fourier spectrum. Section 4 contains a description of the quaternionic phase retrieval problem (QPhR) and the three new algorithms: phase only retrieval, axis only retrieval, and phase axis retrieval. Afterwards, in Sect. 5, we will present the results of the algorithms. Finally, in Sect. 6, we envisage further modifications of QPhR.

Quaternionic setting
A quaternion q is an expression = |q|e uθ = |q|(cos θ + u sin θ) with real-valued a, b, c and d being called components. For i, j and k we have the following multiplication rules: If S(q) = 0, the quaternion is called pure quaternion. The polar form of a quaternion is given in (2). Its components magnitude |q|, axis u and phase θ are determined using the equations [5] The axis u satisfies |u| = 1 and is a unit pure quaternion. A quaternion-valued function f (x, y) can be expressed equivalently to the Cartesian representation of a quaternion (1) as Due to the multiplication of quaternions being non-commutative there are several different versions of QFT in 1D and 2D. Here we apply the 2D single-axis left-sided QFT defined as in [5]: The inverse transformation is defined as Here the so-called axis of QFT µ (not to be confused with the quaternionic axis u) is a unit pure quaternion with μ 2 = −1.

Quaternionic RGB images
RGB images consist of three channels describing the red, green, and blue components of the image. The vector part of a quaternion consists out of three components. Therefore pure quaternions are an ideal way to represent RGB images [5,[8][9][10]. A digital colour image can be described by a discrete function f (n, m), which assigns each pixel [n, m] a Fig. 1 Illustration for the importance of quaternionic Fourier (QF) magnitude, axis and phase: a-c original "Letter" images, d-f reconstructed images, applying d QF magnitude of "G", phase and axis of "R", e QF axis of "G", magnitude and phase of "R", f QF magnitude of "B", axis of "R", phase of "G" quaternion describing the colour at this position: with r , g, b being the red, green, and blue components. The entire RGB image can be represented by a (N × M) matrix of quaternions. The main advantage is that quaternions allow the work with the RGB image as a whole. In consequence, less operations are needed for the QFT of a colour image, since it performs fewer multiplications and additions than three complex Fourier transforms need to transform each colour channel individually; further details can be found in [8,11].

Importance of magnitude, axis and phase
Using Eq. (2) for expressing the Fourier transformed function in polar form we get F(u, v) = |F(u, v)|e v(u,v)ϕ(u,v) with ϕ(u, v) denoting the quaternionic Fourier phase and v(u, v) denoting the quaternionic Fourier axis. Here, the question of the importance of quaternionic Fourier magnitude, axis, and phase spectrum arises. Recapitulating the experiments from Oppenheim et al [12], we adapt them for QFT and colour images in a similar way: (1) Performing QFT for 3 RGB images; (2) In the Fourier domain, the quaternionic Fourier magnitude, axis, and phase are exchanged between the images; (3) Performing inverse QFT for getting RGB images. For exemplification, three images with the individual letters "R", "G" and "B" in red, green, and blue each have been chosen. The experimental results are depicted in Fig. 1. Oppenheim et al [12] noted that the Fourier phase spectrum describes structural information and the fine details of the images. This is also valid for QFT here, even though the object may not have the right colour, as illustrated for the letter "R" in Fig. 1e and "G" in 1f. The QFT magnitude spectrum is relatively negligible in analogy to classical FT [12]. Finally, the QFT axis spectrum v(u, v) as the new part of the quaternionic approach contains the local colour information. In our experiments the quaternionic Fourier axis of the letter "G" is kept unchanged, in consequence, the "G" appears in the right colour in the reconstruction, as shown in Fig. 1e. In the same way, the quaternionic Fourier axis of the letter "R" is maintained in Fig. 1d, f, and results in a correct colour retrieval.

Quaternionic phase retrieval
To describe the QPhR, we follow the description of the classical phase retrieval problem [13][14][15]. A colour image in quaternionic polar form (2) The main difference to the complex case lies in the quaternionic axis u(x, y) as generalisation of the complex unit i. Contrary to i, the axis u consists of three components that are generally unknown and have to be retrieved. Assuming the given magnitudes | f | and |F| in FD, a solution f (x, y) has to fulfil this constraint and prior knowledge, as e.g. positivity of intensity, real-value in the spatial domain for each component or sparsity. It is possible to use projection methods to solve the classical phase retrieval from the view of convex optimization even though the problem itself is not convex [13], or to consider the phase retrieval problem as a minimization problem [16]. This can be applied to the QPhR as well. There are various algorithms, such as the GS algorithm [1], the error reduction algorithm [17], and Fienups algorithms [14,15]. In the following sections we have modified the classical GS algorithm to a quaternionic one, examined with colour images represented by quaternions. All three algorithms follow the scheme in Fig. 2.
They differ in the given components of polar form and the components that will be reconstructed as illustrated by Table 1.

Phase only retrieval
We assume that the magnitudes | f (x, y)| and |F(u, v)| in the spatial domain (SD) and FD are given as well as the axes u(x, y) and v(u, v) in both the object domain and FD. Only the phases θ(x, y) and ϕ(u, v) are to be reconstructed. As starting value for the phase in the FD we use a matrix of random values between −π and π . Together with the known axis v(u, v), we compute e v(u,v)ϕ 1 (u,v) . The differences are the use of the QFT and its inverse instead of the complex FT. Furthermore, not only the magnitudes are replaced with the given values, but also the axes u and v in both domains.

Axis only retrieval
The magnitudes | f (x, y)| and |F(u, v)| in SD and FD are assumed to be known. Furthermore, the phases θ(x, y) and ϕ(u, v) in SD and FD are given. The axes u(x, y) and v(u, v) will be retrieved. Additionally, the prior knowledge of axes being a unit pure quaternion is given. Starting value for the axis is chosen as In contrast to the conventional GS algorithm now the  axes u(x, y) and v(u, v) are normalized, and only their vector parts are used for each iteration.

Phase axis retrieval
Now only the magnitudes | f (x, y)| and |F(u, v)| in SD and FD are given. We want to reconstruct both the axes u(x, y) and v(u, v) and the phases θ(x, y) and ϕ(u, v) in both object and FD. Again we know the axis of a quaternion is a unit pure quaternion. The initial value for the Fourier axis is v 1 = 1 √ 3 (1i + 1 j + 1k) and a constant Fourier phase ϕ 1 (u, v) = π 2 for each quaternion. On both sides now only the resulting magnitudes are replaced with the known magnitudes.

Numerical experiments
All programs were implemented in Matlab R2019a; for quaternionic image representation and the QFT the 'Quaternion and octonion toolbox for Matlab' [18] was applied. The methods were tested on 3 image data sets representing: (1) natural scene images, using CIFAR-10 data, (2) spatially sparse images, using colored MNIST data, (3) textured images, using Brodatz texture data. The results are summarised in Table 2. For illustration, the images "flower", "letters" and "texture" have been chosen here in the following, see Fig. 3. Fig. 3 Example images: a "flower", b "letters" and c "texture" Fig. 4 By phase only retrieval reconstructed images: a "flower", b "letter" and c "texture"

Phase only retrieval
In the phase only retrieval algorithm magnitudes and axes are known and the phase is to be reconstructed. As axis of the QFT µ = 1 √ 3 (1i + 1 j + 1k) was chosen. The algorithm terminates after 200 iterations. Figure 4 shows the results of the phase only retrieval algorithm.
Almost no differences are visible. Some mirror artefacts appear in the image "letters". For an evaluation of the convergence the total square error (SE) for every iteration using a variation of the SE given in [19] is considered  Fig. 5a. The SE of the phases is given by with θ r (n, m) as phase of the reconstructed image and θ o (n, m) as phase of the original image, see Fig. 5b. The quality of the reconstruction regarding visual perception is expressed by the structural similarity index measure (SSIM) between two images X and Y , which is defined [19] as: with μ X and μ Y as values of X , respectively, Y , σ X and σ Y as variances of X and Y , σ XY as covariance of X and Y , c 1 = (k 1 L) 2 and c 2 = (k 2 L) 2 , L is the dynamic range of the  Fig. 10a, E SS I M converges almost to 1, best for the image "flower", reflecting the good reconstruction of natural scene images; for the sparse structure image "letter" and image "texture" E SS I M is only slightly lower.

Axis only retrieval
In the axis only retrieval algorithm magnitudes and phases are assumed to be known and the axes are to be reconstructed. For the QFT the axis µ = 1 √ 2 (1i +1 j +0k) and 200 iteration steps are chosen. The results are shown in Fig. 6. The images appear slightly miscoloured and the background of the image "letters" is also miscoloured since the axis contains information about the colour of an image. In contrast, the general SE (11) converges very fast, as depicted in Fig. 7a. Additionally, mirror artefacts are visible again. The SE for the axes is given as Here, u r (n, m) is the axis of the reconstructed image, while u o (n, m) is the axis of the original image. Again, as shown in Fig. 7b the error converges fast, especially for the sparse structure image "letters" and "texture". Here, E SS I M , which is less than 1, reflects the worse reconstruction, as shown in Fig. 10b.

Phase axis retrieval
Only the magnitudes are given, and both phases and axes have to be reconstructed (applying QFT axis µ = 1 √ 2 (1i +1 j +0k) and 200 iterations). Figure 8 shows the results of the phase axis retrieval algorithm.The images appear greyish, containing only faint colours. The worse reconstruction becomes especially noticeable in the sparse structure image "letters". The background is miscoloured again, the blue colour of the letter "C" is missing, appearing green instead. For all 3 image types the details are retained, indicating a good reconstruction of the phase. The colours appear wrong and mixed up again. Figure 9 shows the quadratic errors (11), (12) and (14). The general SE (Fig. 9a) converges fast again. All SE of the phase (Fig. 9b) converge after few iterations. The SE of the axis (Fig. 9c) decreases, after partly reaching a local maximum first. E SS I M (Fig. 10c) is noticeably below 1, indicating the imperfect reconstruction. The worse result for the "letter" image in axis only retrieval might be caused by the sparse letter structures against white background with colour/mirror artefacts. All methods can reconstruct the phase, while the axis still poses problems since initially not enough constraints for the axis have been included.
6 Further discussions

Influence of the axis
The axis of the Fourier transform µ appears in the kernel of the QFT (8). The only restriction on µ is to be a unit pure quaternion. This means there is an infinite number of possible axes. In the algorithms, we used either µ = 1 The choice of the axis µ is of great importance as it influences the results of algorithms, especially visible in the phase axis retrieval algorithm. Figure 11 shows the results of phase axis retrieval for several choices of µ. Depending on applied µ one of the three primary colours red, green, blue of the letters does not get reconstructed, e.g. in Fig. 11c the originally red "A" appears blue. In Fig. 11f all colours have vanished and the image appears in grayscale. The reason for these effects and the exact influence of the QFT axis µ on the quaternion axis need to be further investigated.

Permutations as constraints
As noted before, not enough constraints are given for a correct reconstruction of the axis of the image. By permuting the colour channels of the image and thus getting two more images artificially, new constraints are created. If the retrieval is executed for all three images at the same time, we can

Extension to other colour spaces and polarimetric images
Experiments in other colour spaces, such as HSV and YCbCr shall be performed. However, the approach must be adapted depending on the structure of the colour space, and the interpretation of quaternionic magnitude, axis and phase therein. Furthermore, our approach might also be applied to other types of (3 × 1)-vector images, e.g. so-called Stokes Vector (SV) images which are measured in polarimetric imaging and can be expressed by their three SV components S 1 , S 2 , S 3 . The magnitudes (expressed at Poincare sphere) constrain S 2 1 + S 2 2 + S 2 3 = S 2 0 , with S 2 0 as total intensity. Assuming local almost continuously increasing accumulated retardations over depth z, then (with r = S 0 , ψ(z) = m 1 z, ξ(z) = S 2 (z) =r cos(ψ) sin(ξ ) = S 0 cos(m 1 z) sin(m 2 z), S 3 (z) =r cos(ψ) cos(ξ ) = S 0 cos(m 1 z) cos(m 2 z).
Now, regarding a birefringent material with locally slowly varying S 1 (x, z) over depth, i.e. for m 1 << m 2 , then S 2 (x, z) can be modelled as the Hilbert transform of S 3 (x, z), which enables an additional prior. Assigning the three SV components to the colour channels yields a pseudo RGB image, as depicted in Fig. 13.

Conclusion and outlook
We have proposed three new QPhR versions, based on the GS algorithm, that reconstruct colour images represented by quaternions. The algorithm retrieves either the phase only, the axis only or both axis and phase. The reconstruction of the axis is most complicated as not enough constraints are given to retrieve the three values that compose the axis. Creating more constraints, e.g. by permuting the colour channels, improves the reconstruction of the algorithms. In the future, the exact convergence rates of the algorithms, especially in relation to each other, should be considered deeper. Additionally, the exact influence of axis µ of the QFT is still to be investigated in greater detail. Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability After requests to authors data and software can be provided.

Conflict of interest No competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.