John von Neumann’s Space-Frequency Orthogonal Transforms

Abstract

Among the invertible orthogonal transforms employed to perform the analysis and synthesis of 2D signals (especially images), the ones defined by means of John von Neumann’s cardinal sinus are extremely interesting. Their definitions rely on transforms similar to those employed to process time-varying 1D signals. This article deals with the extension of John von Neumann’s transforms from 1D to 2D. The approach follows the manner in which the 2D Discrete Fourier Transform was obtained and has the great advantage of preserving the orthogonality property as well as the invertibility. As an important consequence, the numerical procedures to compute the direct and inverse John von Neumann’s 2D transforms can be designed to be efficient thanks to 1D corresponding algorithms. After describing the two numerical procedures, this article focuses on the analysis of their performance after running them on some real-life images. One black and white and one colored image were selected to prove the transforms’ effectiveness. The results show that the 2D John von Neumann’s Transforms are good competitors for other orthogonal transforms in terms of compression intrinsic capacity and image recovery.

1. Introduction

Bidimensional (2D) Signal Processing is an important and large subfield of Signal Processing (SP) [1,2] that encompasses a variety of methods and algorithms. Mathematically, the 2D signals are represented by means or matrices. Typical 2D signals are the still (photographic) images and the complex-valued signals stored in matrices (video streams are rather considered of the 3D or even of multi-dimensional type, where one dimension is always associated with time). Therefore, (still) image processing is a special branch of 2D SP. This article aims to introduce John von Neumann’s Transform (JvNT) devoted to the analysis and synthesis of such signals. The construction of 2D JvNT relies on two other transforms: the (well known) Fourier Transform (FT) and the 1D JvNT, described at length in [3].

Like in the case of 1D signals, the framework of 2D signals constitutes Lebesgue spaces for which both the FT and the scalar product can be defined, as mentioned in [3]. Usually, the bases of 2D signals are generated either by using naturally born 2D waveforms or by coupling two sets of 1D waveforms (one dealing with rows and the other one approaching the columns of a matrix). The orthogonality constraint generally is difficult to verify in the case of 2D waveforms. Therefore, the second approach is often preferred, as the orthogonality of 1D waveforms can easily be preserved.

The extension of the ideal Karhunen–Loève Transform from 1D to 2D is carried out using the first approach, though (see [2]). Instead of working with eigenvalues and eigenvectors of the autocorrelation matrix, like in the case of 1D signals, for 2D signals, some scientists employ singular value decomposition to extract the principal components. Other scientists preserve the 1D approach, but the autocorrelation matrix is computed through multiplication of the signal matrix by the transposed matrix, after extracting the average.

In the case of transforms based on orthogonal polynomials (as mentioned in [3]), the second approach is adopted to jump into the 2D signal spaces. Nevertheless, the complexity of such transforms can tremendously increase, as shown, for example, in [4] and, very recently, in [5].

Beside the ideal and the polynomial-based transforms, the world of 2D orthogonal transforms can be structured into four classes, like in the case of 1D transforms [3], by replacing the time domain with the space domain: harmonic, space-frequency, space-scale, space-frequency-scale. Basically, any 1D transform can be extended to a 2D transform, under the same name.

The harmonic class is naturally generated by the FT. The extension from 1D to 2D is performed using the “two 1D sets” approach, which will be shortly reviewed within the next section. Despite its age, this transform continues to be employed in applications (although some approaches reported nowadays in the literature are rather naïve). For example, in [6], the 2D FT symmetry properties are analyzed. Symmetry is an important and desirable property of transforms, and can help reduce the computational burden. (Unfortunately, not all transforms are symmetric or conjugate-symmetric like FT. Fortunately, both 1D and 2D JvNTs are verified to possess such properties.) Besides 2D FT, two other versions have been defined: Windowed Fourier Transform (WFT, which is quite old, too) [7] and fractional Fourier Transform (frFT, which was quite recently introduced) [8,9]. In many applications, and especially in modern signal compression technology, 2D FT was replaced by the family of 2D Discrete Cosine Transforms, for which the analysis coefficients are real-valued (and not complex-valued, like in the case of FT). Note, however, that the cosine coefficients are not necessarily symmetric. Recent publications, like [10,11,12,13], prove that the employment of these transforms still is quite intense.

Space-frequency transforms also are obtained by windowing the signals. Sometimes, windowed transforms are referred to as smooth wavelet transforms. JvNTs are actually such transforms, where the working window is the John von Neumann (JvN) cardinal sinus. Besides JvNTs, some other remarkable transforms are reported in the literature. For example, the old Morlet wavelet (well known as the Mexican hat), which initially was devoted to geological prospections, is employed nowadays in image compression [14]. Another example is given by the Gauss–Gabor wavelet, which was extended to 2D signals in [15] using an interesting approach.

The class of space-scale transforms is quite large and, usually, is split into two groups: empirical transforms and generic transforms. Fractal 2D signals are mainly processed with such transforms. In fact, the fractal wavelets constitute the foundation of this class. The old Haar wavelet (also known as the lazy wavelet), as a typical empirical transform, works very well with images, as recently reported in [16,17]. The more sophisticated Walsh–Hadamard Transform (still empirical, though) can also be adapted to process images—see, for example, [18,19] (where it works in combination with Haar wavelet). Another empirical transform with increasing impact in image processing is the (less known) Slant Transform [20]. Although the extension of this transform to 2D signals is not easy, it was employed to denoise images, like in [21], and to process video streams (among other transforms), as recently reported in [22].

Without any doubt, nowadays, generic (fractal) orthogonal and biorthogonal wavelets are the most employed tools in image processing. Since their inception by Meyer, Mallat and Daubechies (more than 35 years ago), the literature reporting how these wavelets work in conjunction with images has become so vast that any attempt to encompass all the sound (and sometimes amazing) results is very likely doomed to fail. One can only cite a few interesting recent works, such as [23] (with application to image filtering), [24] (in which principal component analysis is developed by means of wavelets, aiming to achieve image denoising), [25] (where image fusion is performed) and [26] (which deals with the efficient implementation of 2D wavelet transforms).

Transforms of the space-frequency-scale class exhibit great complexity and usually are non-orthogonal. Basically, any space-frequency or space-scale transform can evolve towards this class (including JvNT). The basis of the Lebesgue space is generated by means of all three operators: time-shifting, frequency modulation and scaling. They are applied (in this order) to a mother waveform (or wavelet) and generate a whole dictionary of space-frequency-scale atoms. Analysis of orthogonality and invertibility of the transform based on such a dictionary is a difficult task (and an open problem). Working with non-orthogonal bases is, however, necessary in many applications where image interpretation is required (like in medicine, geography, and underwater or space exploration). In this case, in opposition with image compression or denoising, keeping a high level of image redundancy is crucial.

By determining the difference from 1D transforms, the collection of 2D transforms was enriched with hybrid transforms, born through a combination of two or more already known transforms. For example, one can find a joint Cosine and Karhunen–Loève Transform in [27], a joint Cosine and Hadamard Transform in [28], and a joint Walsh and fractional Fourier Transform in [29].

The exploration of the scientific literature performed within this introduction is far from complete. But the goal was not to review the entire 2D SP subfield. This study aimed to correctly place JvNT in the wide world of orthogonal transforms and to show that, to the best of our knowledge, no other article in the literature has reported a similar approach and similar results to this article.

The following sections are presented next. In Section 2, after a short overview of 1D JvNT (already introduced in [3]), the extension to space signals is described. We prove that working with two 1D dictionaries is equivalent to working with a naturally born 2D dictionary. Section 3 is devoted to the design of the analysis and synthesis of numerical algorithms that allow for the efficient implementation of 2D JvNT. The simulation results and a discussion of the obtained results can be found in Section 4. The algorithms of Section 3 were implemented and tested on two types of images: black and white, and color. Some concluding remarks, an acronym list and a reference list complete the article.

2. Theoretical Background

2.1. Short Overview of JvNT for Discrete Time 1D Signals

The use of JvNT for 1D signals was described at length in [3]. Within this section, only a short overview of main notations and results from [3] is presented. The starting point is the cardinal sinus below (which the great scientist JvN defined long time ago [30]):

$$\mathsf{\nu }(t)=\mathrm{sinc}(\mathsf{\pi }t)=\frac{\sin (\mathsf{\pi }t)}{\mathsf{\pi }t}, \forall t\in \mathbb{R}.$$

(1)

A dictionary of time-frequency atoms (tfas) can be built by applying time-shifting, frequency modulation and discretization to Function (1). More specifically, the generic tfas can be expressed as follows:

$${\mathsf{\nu }}_{T_s}^{[p,k]}[n]={\mathrm{e}}^{2\mathsf{\pi }kn{T}_s\mathrm{j}}\mathrm{sinc}\left[\mathsf{\pi }\left(n{T}_s-p\right)\right], \forall n,p\in \mathbb{Z}, \forall k\in \mathbb{N},$$

(2)

where $T_s\in (0,1]$ is the sampling period, $n$ is the discrete (normalized) time, $p$ is the time-shifting index and $k$ is the frequency modulation index. In [3], it is proven that if $T_s=1/K$, where $K\in {\mathbb{N}}^{\ast }$ is a sampling frequency (with integer values measured in Hz), then the reduced dictionary ${\mathsf{V}}_K={\left\{{\mathsf{\nu }}_{T_s}^{[p,k]}\right\}}_{p\in \mathbb{Z},k\in \overline{0,K-1}}$ is an orthogonal basis of the Lebesgue–Hilbert space ${\mathsf{l}}^2$ (of discrete time signals). In this case, the notation ${\mathsf{\nu }}_{T_s}^{[p,k]}$ changes to ${\mathsf{\nu }}_K^{[p,k]}$, and Equation (2) becomes

$${\mathsf{\nu }}_K^{[p,k]}[n]={\mathrm{e}}^{\frac{2nk\mathsf{\pi }}{K}\mathrm{j}}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{n}{K}-p\right)\right], \forall n,p\in \mathbb{Z}, \forall k\in \overline{0,K-1}.$$

(3)

Consequently, the discrete JvNT is obtained by using the orthogonal basis ${\mathsf{V}}_K$. Thus, the JvN coefficients of the signal $x\in {\mathsf{l}}^2$ can be computed by projecting it on each tfa:

$$\begin{gathered} {\mathsf{N}}_K(x)[p,k]=X_K[p,k]=〈x,{\mathsf{\nu }}_K^{[p,k]}〉={\displaystyle \sum _{n\in \mathbb{Z}}x[n]\overline{{\mathsf{\nu }}_K^{[p,k]}[n]}}={\displaystyle \sum _{n\in \mathbb{Z}}x[n]\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{n}{K}-p\right)\right]{\mathrm{e}}^{-\frac{2kn\mathsf{\pi }}{K}\mathrm{j}}}, \\ \forall p\in \mathbb{Z}, \forall k\in \overline{0,K-1}. \end{gathered}$$

(4)

To recover the signal, expansion on the basis of ${\mathsf{V}}_K$ with coefficient (4) is employed:

$$x[n]=\frac{1}{K}{\displaystyle \sum _{p\in \mathbb{Z}}{\displaystyle \sum _{k=0}^{K-1}{X}_K[p,k]{\mathsf{\nu }}_K^{[p,k]}[n]}}=\frac{1}{K}{\displaystyle \sum _{p\in \mathbb{Z}}{\displaystyle \sum _{k=0}^{K-1}{X}_K[p,k]\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{n}{K}-p\right)\right]{\mathrm{e}}^{\frac{2kn\mathsf{\pi }}{K}\mathrm{j}}}}, \forall n\in \mathbb{Z}.$$

(5)

If the signal $x\in {\mathsf{l}}^2$ is real-valued, the following remarkable property is verified:

$$X_K[p,K-k]=\overline{X_K[p,k]}, \forall p\in \mathbb{Z}, \forall k\in \overline{0,K-1}.$$

(6)

This means the JvN coefficients are conjugate-symmetric, which allows for a reduction in the computational burden by only evaluating approximately half of them. Moreover, the synthesis Equation (5) can benefit from symmetry in implementation.

In Definition (4) and Equation (5), signals from ${\mathsf{l}}^2$ either have infinite or finite support. In the second case, signals not only have finite energy, but are also stable (absolutely summable), i.e., they also belong to ${\mathsf{l}}^1\subset {\mathsf{l}}^2$. Consequently, Fourier Transform (FT), and especially Discrete Fourier Transform (DFT), can be computed for such signals.

Also, for finite support signals, the JvNT yields when working with a finite number of JvN coefficients, as shown in [3]. More specifically, depending on the signal length $N$ and accuracy threshold $\mathsf{\epsilon }>0$ (which allows for truncation of the signal (1) at the user’s convenience), the lower and upper bounds of time-shifting index $p$ can be derived: $p\in \overline{P_{\min },P_{\max }}$. The upper limit of the frequency modulation index, $K$, depends on $N$, as well, and can be set according to the Gabor–Heisenberg Uncertainty Principle [2], such that the time and frequency resolutions of the JvN spectrogram be approximately equal. For example, this requirement is fulfilled if $K=⌊\sqrt{2N}+0.5⌋$.

Although the synthesis of a finite-length signal with truncated tfas is not exact (except for the normalized instants located at integer multiples of $K$), the reconstruction accuracy can be controlled by the user through the parameters $\mathsf{\epsilon }$ and $K$.

The extension of JvNT to spatial signals (images) strongly relies on the properties above of JvNT devoted to 1D signals.

2.2. JvNT for Discrete-Space 2D Signals

Although it is possible to develop a theory of continuous-space JvNT for 2D signals, the practical importance of such a transform is rather small. Nowadays, 2D signals (especially images) are digitally represented, processed and transmitted. Therefore, the following development only concerns digital signals, seen as (still) images.

Let $A\in {\mathbb{R}}^{NR\times NC}$ be a matrix playing the role of a 2D discrete spatial signal. The matrix has $NR$ rows and $NC$ columns. If $A$ is a digital representation of some (still) black and white (bw) image, then its elements are referred to as pixels. Hereafter, by convention, the term of pixel will be associated with the current element $a_{nr,nc}$ of matrix $A$, located at coordinates $(nr,nc)$ (for $nr\in \overline{1,NR}$ and $nc\in \overline{1,NC}$). In SP, such a matrix can be transformed in two ways: by means of a dictionary with 2D genuinely born atoms or by means of two dictionaries with 1D atoms.

In the first approach, usually, the atoms of the dictionary are obtained by spatially rotating a 1D mother waveform (mw). For example, in the case of JvN mw (1) and the matrix $A$, one (non-unique) rotated version is as follows:

$$\mathsf{\nu }(x,y)=\mathrm{sinc}\left(\mathsf{\pi }\sqrt{\frac{2NR}{NR+NC}{x}^2+\frac{2NC}{NR+NC}{y}^2}\right), \forall x,y\in \mathbb{R}.$$

(7)

If $A$ is a square matrix ($NR=NC$), then the mw (7) becomes

$$\mathsf{\nu }(x,y)=\mathrm{sinc}\left(\mathsf{\pi }\sqrt{x^2+y^2}\right), \forall x,y\in \mathbb{R}.$$

(8)

In the top image in Figure 1, the shape of mw (8) (for $NR=NC$) is drawn, while at the bottom, 2D FT is represented (spectrogram and phase surface).

The imperfections in the displayed spectrogram are due to the truncation of mw. In fact, the two dominant corners of the spectral surface are cylindrical. This property sensibly complicates the analysis of orthogonality. Such an analysis, although interesting, is not developed within this article.

The second approach is more promising due to the strong similarity with 2D FT. Recall that the 2D DFT of matrix $A$ is defined as below [1,2]:

$$\mathsf{F}(A)\left({\mathsf{\omega }}_x,{\mathsf{\omega }}_y\right)={\displaystyle \sum _{nr=0}^{NR-1}{\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}{\mathrm{e}}^{-\mathrm{j}\left(nr{\mathsf{\omega }}_x+nc{\mathsf{\omega }}_y\right)}}}={\displaystyle \sum _{nr=0}^{NR-1}{\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}{\mathrm{e}}^{-\mathrm{j}nr{\mathsf{\omega }}_x}{\mathrm{e}}^{-\mathrm{j}nc{\mathsf{\omega }}_y}}}, \forall {\mathsf{\omega }}_x,{\mathsf{\omega }}_y\in [-\mathsf{\pi },+\mathsf{\pi }].$$

(9)

The last expression in definition (9) allows us to write the following matrix form:

$$\mathsf{F}(A)\left({\mathsf{\omega }}_x,{\mathsf{\omega }}_y\right)=\overline{e_{NR}^T\left({\mathsf{\omega }}_x\right)}A\overline{e_{NC}\left({\mathsf{\omega }}_y\right)}, \forall {\mathsf{\omega }}_x,{\mathsf{\omega }}_y\in [-\mathsf{\pi },+\mathsf{\pi }],$$

(10)

where, by definition:

$$e_N^T\left(\mathsf{\omega }\right)=\left[1\hspace{1em} {\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\hspace{1em} {\mathrm{e}}^{2\mathrm{j}\mathsf{\omega }}\hspace{1em} \dots \hspace{1em} {\mathrm{e}}^{\mathrm{j}(N-1)\mathsf{\omega }}\right], \forall \mathsf{\omega }\in [-\mathsf{\pi },+\mathsf{\pi }].$$

(11)

Equation (10) reveals that two harmonic vectors can be employed to compute the FT: one operating on matrix columns and another one operating on matrix rows. This property suggests that two dictionaries of orthogonal 1D tfas can be built (3): one associated with matrix columns and another one associated with matrix rows. Each dictionary is configured by some sampling rate, namely $KR$ and $KC$, respectively. It is not necessary that the parameters $KR$ and/or $KC$ be directly correlated with the matrix sizes $NR$ and $NC$. Nevertheless, according to the previous subsection, one can set $KR=⌊\sqrt{2NR}+0.5⌋$ and $KC=⌊\sqrt{2NC}+0.5⌋$. Also, the time-shifting index $p$ of the previous section becomes, in this context, a space-shifting index. Since the matrix has finite sizes, the bounds of the space-shifting indices can be evaluated for rows, as well as for columns, and they are not necessarily the same. They are denoted as $P{R}_{\min }$ and $P{R}_{\max }$ for columns (of length $NR$) and by $P{C}_{\min }$ and $P{C}_{\max }$ for rows (of length $NC$). Usually, the accuracy threshold $\mathsf{\epsilon }$ is unique (for example, set to 1%).

Definition 1.

The practical JvNT of 2D signals is defined as follows:

$$\begin{gathered} {\mathsf{N}}_{KR,KC}(A)[pr,kr,pc,kc]={\left(\overline{{\mathsf{\nu }}_{KR}^{[pr,kr]}}\right)}^{T}A\overline{{\mathsf{\nu }}_{KC}^{[pc,kc]}}, \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}, \end{gathered}$$

(12)

where

$$\begin{gathered} \left\{\begin{array}{l}{\mathsf{\nu }}_{KR}^{[pr,kr]}={\left[\begin{array}{c}{\mathsf{\nu }}_{KR}^{[pr,kr]}[0]\hspace{1em}{\mathsf{\nu }}_{KR}^{[pr,kr]}[1]\hspace{1em}{\mathsf{\nu }}_{KR}^{[pr,kr]}[2]\hspace{1em}\cdots \hspace{1em}{\mathsf{\nu }}_{KR}^{[pr,kr]}[NR-1]\end{array}\right]}^T\\ {\mathsf{\nu }}_{KC}^{[pc,kc]}={\left[\begin{array}{c}{\mathsf{\nu }}_{KC}^{[pc,kc]}[0]\hspace{1em}{\mathsf{\nu }}_{KC}^{[pc,kc]}[1]\hspace{1em}{\mathsf{\nu }}_{KC}^{[pc,kc]}[2]\hspace{1em}\cdots \hspace{1em}{\mathsf{\nu }}_{KC}^{[pc,kc]}[NC-1]\end{array}\right]}^T\end{array}\right., \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(13)

Although not obvious, the definition above is equivalent to working in 2D signal space with a 2D JvN function:

$${\mathsf{\nu }}_{2\mathrm{D}}(x,y)=\mathrm{sinc}(\mathsf{\pi }x)\mathrm{sinc}(\mathsf{\pi }y)=\mathsf{\nu }(x)\mathsf{\nu }(y), \forall x,y\in \mathbb{R}.$$

(14)

Recall that Function (1) is also the impulse response of an ideal Low-Pass Filter (LPF). This property is inherited by the discrete signal ${\mathsf{\nu }}_{T_s}^{[0,0]}$, regardless of the value of the sampling period $T_s\in (0,1]$. More specifically, assume the following ideal digital LPF has a frequency response expressed as

$$H_1\left({\mathrm{e}}^{\mathrm{j}\mathsf{\omega }}\right)=\left\{\begin{array}{l}1\hspace{1em},\mathsf{\omega }\in \left[-{\mathsf{\omega }}_c,{\mathsf{\omega }}_c\right]\hfill \\ 0\hspace{1em},\mathsf{\omega }\in [-\mathsf{\pi },+\mathsf{\pi }]\backslash \left[-{\mathsf{\omega }}_c,{\mathsf{\omega }}_c\right]\hfill \end{array}\right.={\chi }_{\left[-{\mathsf{\omega }}_c,{\mathsf{\omega }}_c\right]}(\mathsf{\omega }), \forall \mathsf{\omega }\in [-\mathsf{\pi },+\mathsf{\pi }],$$

(15)

where ${\chi }_{\left[a,b\right]}$ is the index function of interval $\left[a,b\right]$ and ${\mathsf{\omega }}_c=\mathsf{\pi }{T}_s\in (0,\mathsf{\pi })$ is the cut-off normalized pulsation. The impulse response of the filter is then

$$h_1[n]=\frac{{\mathsf{\omega }}_c}{\mathsf{\pi }}{\mathsf{\nu }}_{T_s}^{[0,0]}\left(\frac{{\mathsf{\omega }}_c}{\mathsf{\pi }}n\right), \forall n\in \mathbb{Z}.$$

(16)

In the 2D space of discrete signals, Equations (15) and (16) become

$$\begin{array}{l}{H}_2\left({\mathrm{e}}^{\mathrm{j}{\mathsf{\omega }}_x},{\mathrm{e}}^{\mathrm{j}{\mathsf{\omega }}_y}\right)=\left\{\begin{array}{l}1\hspace{1em},\left({\mathsf{\omega }}_x,{\mathsf{\omega }}_y\right)\in \left[-{\mathsf{\omega }}_{cx},{\mathsf{\omega }}_{cx}\right]\times \left[-{\mathsf{\omega }}_{cy},{\mathsf{\omega }}_{cy}\right]\hfill \\ 0\hspace{1em},\mathrm{otherwise}\hfill \end{array}\right.={\chi }_{\left[-{\mathsf{\omega }}_{cx},{\mathsf{\omega }}_{cx}\right]}\left({\mathsf{\omega }}_x\right){\chi }_{\left[-{\mathsf{\omega }}_{cy},{\mathsf{\omega }}_{cy}\right]}\left({\mathsf{\omega }}_y\right)=\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =H_1\left({\mathrm{e}}^{\mathrm{j}{\mathsf{\omega }}_x}\right)H_1\left({\mathrm{e}}^{\mathrm{j}{\mathsf{\omega }}_y}\right),\forall {\mathsf{\omega }}_x,{\mathsf{\omega }}_y\in [-\mathsf{\pi },+\mathsf{\pi }];\end{array}$$

(17)

$$h_2[nx,ny]=\frac{{\mathsf{\omega }}_{cx}{\mathsf{\omega }}_{cy}}{{\mathsf{\pi }}^2}{\mathsf{\nu }}_{2\mathrm{D}}\left(\frac{{\mathsf{\omega }}_{cx}}{\mathsf{\pi }}nx,\frac{{\mathsf{\omega }}_{cy}}{\mathsf{\pi }}ny\right)=\frac{{\mathsf{\omega }}_{cx}{\mathsf{\omega }}_{cy}}{{\mathsf{\pi }}^2}\mathsf{\nu }\left(\frac{{\mathsf{\omega }}_{cx}}{\mathsf{\pi }}nx\right)\mathsf{\nu }\left(\frac{{\mathsf{\omega }}_{cy}}{\mathsf{\pi }}ny\right), \forall nx,ny\in \mathbb{Z}.$$

(18)

In Equations (17) and (18), ${\mathsf{\omega }}_{cx}=\mathsf{\pi }{X}_s$ and ${\mathsf{\omega }}_{cy}=\mathsf{\pi }{Y}_s$, where $X_s\in (0,1]$ and $Y_s\in (0,1]$ are sampling distances along the Ox and Oy axes, respectively. Thus, in fact, Equation (18) is equivalent to:

$$h_2[nx,ny]=\frac{{\mathsf{\omega }}_{cx}{\mathsf{\omega }}_{cy}}{{\mathsf{\pi }}^2}{\mathsf{\nu }}_{X_s}^{[0,0]}\left(\frac{{\mathsf{\omega }}_{cx}}{\mathsf{\pi }}n\right){\mathsf{\nu }}_{Y_s}^{[0,0]}\left(\frac{{\mathsf{\omega }}_{cy}}{\mathsf{\pi }}n\right)=h_1[nx]h_1[ny], \forall nx,ny\in \mathbb{Z}.$$

(19)

While in the space of 1D signals, the LPF frequency characteristic is a rectangle over the $[-{\mathsf{\omega }}_c,{\mathsf{\omega }}_c]$ band, and in the space of 2D signals, it is a parallelepiped over the $\left[-{\mathsf{\omega }}_{cx},{\mathsf{\omega }}_{cx}\right]\times \left[-{\mathsf{\omega }}_{cy},{\mathsf{\omega }}_{cy}\right]$ band (rectangular). Orthogonality can be analyzed by shifting the parallelepiped along the pulsation axes such that the rectangular bands are almost disjunct (i.e., they intersect most along a vertical plane or a vertical line).

The only problem is that of defining appropriate scalar products in discrete space and frequency. In general, defining scalar products between 2D signals (i.e., between matrices) is a difficult task. Fortunately, as Equations (17) and (19) clearly reveal, both the impulse response and the frequency response can be expressed in factorial form with the help of 1D signals, one for each space or frequency axes. Hence, one can write

$$h_2\equiv h_1^{\mathrm{Ox}}{h}_1^{\mathrm{Oy}}, H_2\equiv H_1^{\mathrm{Ox}}{H}_1^{\mathrm{Oy}},$$

(20)

with natural notations. Identity (20) allows us to define special forms of scalar products in the context of 2D signals with the help of scalar products already defined in the context of 1D signals:

$$\left\{\begin{array}{l}〈h_2,g_2〉=〈h_1^{\mathrm{Ox}}{h}_1^{\mathrm{Oy}},g_1^{\mathrm{Ox}}{g}_1^{\mathrm{Oy}}〉=〈h_1^{\mathrm{Ox}},g_1^{\mathrm{Ox}}〉〈h_1^{\mathrm{Oy}},g_1^{\mathrm{Oy}}〉\\ 〈H_2,G_2〉=〈H_1^{\mathrm{Ox}}{H}_1^{\mathrm{Oy}},G_1^{\mathrm{Ox}}{G}_1^{\mathrm{Oy}}〉=〈H_1^{\mathrm{Ox}},G_1^{\mathrm{Ox}}〉〈H_1^{\mathrm{Oy}},G_1^{\mathrm{Oy}}〉\end{array}\right..$$

(21)

Definition (21), together with Definition (17) and Equation (19), enable the analysis of the orthogonality of 2D tfas by means of the orthogonality of corresponding 1D tfas. All the orthogonality and invertibility results proven in [3] for the 1D time-frequency dictionaries of JvN atoms can straightforwardly be transferred to 2D space-frequency dictionaries. The generic space-frequency atom (sfa) of the JvN 2D dictionary is as follows:

$$\begin{gathered} {\mathsf{\nu }}_{KR,KC}^{[pr,pc,kr,kc]}[nr,nc]={\mathsf{\nu }}_{KR}^{[pr,kr]}[nr]{\mathsf{\nu }}_{KC}^{[pr,kr]}[nc], \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(22)

Thanks to the factorization exhibited in Equation (22), the projection of matrix $A$ on any sfa of the dictionary can be defined as follows:

$$\begin{gathered} 〈A,{\mathsf{\nu }}_{KR,KC}^{[pr,pc,kr,kc]}[nr,nc]〉={\displaystyle \sum _{nr=0}^{NR-1}{\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}{\mathsf{\nu }}_{KR}^{[pr,kr]}[nr]{\mathsf{\nu }}_{KC}^{[pr,kr]}[nc]}}, \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(23)

At the same time, Definition (12) can be expressed at length:

$$\begin{gathered} {\mathsf{N}}_{KR,KC}(A)[pr,kr,pc,kc]={\displaystyle \sum _{nr=0}^{NR-1}{\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}{\mathsf{\nu }}_{KR}^{[pr,kr]}[nr]{\mathsf{\nu }}_{KC}^{[pr,kr]}[nc]}}=〈A,{\mathsf{\nu }}_{KR,KC}^{[pr,pc,kr,kc]}[nr,nc]〉, \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(24)

From Equation (24), it results that the analysis of 2D signal $A$ is equivalently performed by using a couple of 1D dictionaries (according to Definition (12)) or with the help of sfas taken from a 2D dictionary (according to Definition (23)).

Furthermore, Equation (24) can be expressed in implementation form below:

$$\begin{gathered} {\mathsf{N}}_{KR,KC}(A)[pr,kr,pc,kc]={\displaystyle \sum _{nr=0}^{NR-1}{\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nr}{KR}-pr\right)\right]}}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nc}{KC}-pc\right)\right]{\mathrm{e}}^{-2\mathsf{\pi }\left(\frac{krnr}{KR}+\frac{kcnc}{KC}\right)\mathrm{j}}, \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(25)

The JvN coefficients (25) can be grouped into a 4D array and cannot be represented visually. To compute them, the previously mentioned symmetry property can be exploited for each harmonic index, provided that the 2D signal is real-valued.

As already mentioned, the orthogonality of the 2D dictionary involves the orthogonality of the transform (12). Consequently, the practical inverse of JvNT can be computed like in [3]. More specifically, the extension to spatial signals is expressed as follows:

• if $nr=mr·KR\in KR\mathbb{N}\cap \overline{0,NR-1}$ and $nc=mc·KC\in KC\mathbb{N}\cap \overline{0,NC-1}$:

  $$\tilde{a}[mrKR+1,mcKC+1]=\frac{1}{KR·KC}{\displaystyle \sum _{kr=0}^{KR-1}{\displaystyle \sum _{kc=0}^{KC-1}{\mathsf{N}}_{KR,KC}(A)[mr,kr,mc,kc]}};$$

  (26)
• if $nr=mr·KR\in KR\mathbb{N}\cap \overline{0,NR-1}$ and $nc\in \overline{0,NC-1}\backslash KC\mathbb{N}$:

  $$\tilde{a}[mrKR+1,nc+1]=\frac{{\displaystyle \sum _{pc=P{C}_{\min }}^{P{C}_{\max }}{\displaystyle \sum _{kc=0}^{KC-1}\left({\displaystyle \sum _{kr=0}^{KR-1}{\mathsf{N}}_{KR,KC}(A)[mr,kr,pc,kc]}\right){\mathsf{\nu }}_{KC}^{[pc,kc]}[nc]}}}{KR·KC{\displaystyle \sum _{pc=P{C}_{\min }}^{P{C}_{\max }}{\left({\mathsf{\nu }}_{KC}^{[pc,0]}[nc]\right)}^2}};$$

  (27)
• if $nr\in \overline{0,NR-1}\backslash KR\mathbb{N}$ and $nc=mc·KC\in KC\mathbb{N}\cap \overline{0,NC-1}$:

  $$\tilde{a}[nr+1,mcKC+1]=\frac{{\displaystyle \sum _{pr=P{R}_{\min }}^{P{R}_{\max }}{\displaystyle \sum _{kr=0}^{KR-1}\left({\displaystyle \sum _{kc=0}^{KC-1}{\mathsf{N}}_{KR,KC}(A)[pr,kr,mc,kc]}\right){\mathsf{\nu }}_{KC}^{[pc,kc]}[nc]}}}{KR·KC{\displaystyle \sum _{pr=P{R}_{\min }}^{P{R}_{\max }}{\left({\mathsf{\nu }}_{KR}^{[pr,0]}[nr]\right)}^2}};$$

  (28)
• if $nr\in \overline{0,NR-1}\backslash KR\mathbb{N}$ and $nc\in \overline{0,NC-1}\backslash KC\mathbb{N}$:

  $$\tilde{a}[nr+1,nc+1]=\frac{{\displaystyle \sum _{pr=P{R}_{\min }}^{P{R}_{\max }}{\displaystyle \sum _{kr=0}^{KR-1}{\displaystyle \sum _{pc=P{C}_{\min }}^{P{C}_{\max }}{\displaystyle \sum _{kc=0}^{KC-1}{\mathsf{N}}_{KR,KC}(A)[pr,kr,pc,kc]{\mathsf{\nu }}_{KR}^{[pr,kr]}[nr]{\mathsf{\nu }}_{KC}^{[pc,kc]}[nc]}}}}}{KR·KC\left[{\displaystyle \sum _{pr=P{R}_{\min }}^{P{R}_{\max }}{\left({\mathsf{\nu }}_{KR}^{[pr,0]}[nr]\right)}^2}\right]\left[{\displaystyle \sum _{pc=P{C}_{\min }}^{P{C}_{\max }}{\left({\mathsf{\nu }}_{KC}^{[pc,0]}[nc]\right)}^2}\right]}.$$

  (29)

(The notations $KR\mathbb{N}$ and $KC\mathbb{N}$ stand for the sets of all non-negative integer multiples of numbers $KR$ and $KC$, respectively).

The reconstructed matrix $\tilde{A}={\left[{\tilde{a}}_{nr,nc}\right]}_{\begin{array}{l}nr\in \overline{1,NR}\\ nc\in \overline{1,NC}\end{array}}$ approximates the original matrix $A$. Since the threshold $\mathsf{\epsilon }$ is constant, the accuracy of approximation $\tilde{\mathbf{A}}$ only depends on the two sampling rates $KR$ and $KC$, which are set according to the Uncertainty Principle. To evaluate the accuracy, the following cost function can be employed:

$$\mathsf{A}(KR,KC,\mathsf{\epsilon })=\frac{100}{1+10\frac{{‖A-\tilde{A}‖}_2}{{‖A‖}_2}} [\%], \forall KR,KC\in {\mathbb{N}}^{\ast }, \forall \mathsf{\epsilon }>0,$$

(30)

where ${‖A‖}_2$ is the 2-norm of matrix $A$ (the largest singular value). Definition (30) employs the function $1/(1+x)$, which compresses the interval $[0,+\infty )$ into the interval $(0,1]$. The cost Function (30) is often referred to as fitness.

Equations (26)–(29) can involve significant computational burden, depending on the size of the original matrix. To reduce the runtime of the corresponding numerical algorithm, the symmetry properties should be accounted for at the expense of higher programming effort. Working with arrays with more than two dimensions is time-consuming within almost all programing environments. Therefore, using any property that can reduce the computational burden is worthwhile.

3. Numerical Algorithms to Implement 2D JvNT

Two algorithms were designed and implemented based on the previous section. Only the case of real-valued matrices was considered for direct 2D JvNT, because complex-valued matrices are unusual in real-life applications. Nevertheless, an algorithm for complex-valued matrices can be designed according to general definition (4). However, another procedure was designed to implement the inverse 2D JvNT. Although the genuine matrix to achieve approximately recovery is real-valued, the 4D array of JvN coefficients includes complex-valued values. In this case, Equations (26)–(29) were employed to design the corresponding algorithm.

Both algorithms were designed to allow implementation within the MATLAB™ programming environment (version 2018), although other high-level languages can be employed, as well (such as C++ or Python).

3.1. Direct JvNT Algorithm for Real-Valued 2D Signals

The analysis algorithm requires some preliminary preparations. The aim is to design a procedure that invokes Algorithm 1 from [3]. In this way, most of the computations are performed just like in the case of 1D signals. Thus, the symmetry property can fully be exploited (the direct implementation of Equation (4) is also possible, even when considering the symmetry).

The computations in Equation (4) can be organized by using the matrix formalism. Thus, the equivalent expression of JvN analysis coefficients is as follows:

$$\begin{gathered} {\mathsf{N}}_{KR,KC}(A)[pr,kr,pc,kc]=\underset{\mathrm{TvN}1\mathrm{D} -ℂ}{\underbrace{{\displaystyle \sum _{nr=0}^{NR-1}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nr}{KR}-pr\right)\right]{\mathrm{e}}^{-\frac{2\mathsf{\pi }krnr}{KR}\mathrm{j}}\underset{\mathrm{TvN}1\mathrm{D} -\mathbb{R}}{\underbrace{\left({\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nc}{KC}-pc\right)\right]{\mathrm{e}}^{-\frac{2\mathsf{\pi }kcnc}{KC}\mathrm{j}}}\right)}}}}} \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(31)

In Equation (31), two JvN 1D Transforms are emphasized: one that operates with real-valued signals and another one that deals with complex-valued signals. Focus first on the inner JvNT. Practically every row of the real-valued matrix $A$ produces a matrix with complex-valued coefficients:

$$\begin{gathered} c_{KC}^{nr+1}[pc-P{C}_{\min }+1,kc+1]={\mathsf{N}}_{KC}\left({\mathrm{row}}_{nr+1}(A)\right)[pc,kc]={\displaystyle \sum _{nc=0}^{NC-1}{a}_{nr+1,nc+1}\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nc}{KC}-pc\right)\right]{\mathrm{e}}^{-\frac{2\mathsf{\pi }kcnc}{KC}\mathrm{j}}}, \\ \forall nr\in \overline{0,NR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(32)

All matrices from Equation (32), ${\left\{C_{KC}^{nr+1}\right\}}_{nr\in \overline{0,NR-1}}$, can be packed into a 3D array, $C_{KC}$, where each matrix is a layer. The sizes of the array $C_{KC}$ are as follows: $PC=P{C}_{\max }-P{C}_{\min }+1$ rows, $KC$ columns and $NR$ layers. Unfortunately, the second JvNT does not act on the columns of each layer, but on vectors extracted from the array $C_{KC}$ along the layers. This is because the outer JvNT in Equation (31) is configured to work with $NR$-length signals. Thus, the 3D array should be reshaped into a large matrix, say $D_{KC}$, with $NR$ rows and $M=PC·KC$ columns. The matrix $D_{KC}$ can be obtained by linearizing each matrix $C_{KC}^{nr+1}$ along its rows and stacking the resulting long rows from top to bottom. Linearization can be implemented by using the well-known Theorem of Division with Remainder (TDR) between integers. Hence, the current element of the matrix $D_{KC}$ is:

$$d_{KC}[nr+1,m+1]=c_{KC}^{nr+1}\left[⌊\frac{m}{KC}⌋+1,m\%KC+1\right], \forall nr\in \overline{0,NR-1}, \forall m\in \overline{0,M-1},$$

(33)

where $⌊a⌋$ is the (lower) integer part of number $a\in \mathbb{R}$, while $n\%N$ stands for the reminder of division $n/N$ (with $n\in \mathbb{Z}$ and $N\in {\mathbb{N}}^{\ast }$).

Now, the 1D JvNT can be applied to each column of the complex-valued matrix $D_{KC}$ (for which the symmetry property cannot be exploited):

$$\begin{gathered} c_{KR,KC}^{m+1}[pr-P{R}_{\min }+1,kr+1]={\mathsf{N}}_{KR}\left({\mathrm{col}}_{m+1}\left(D_{KC}\right)\right)[pr,kr]={\displaystyle \sum _{nr=0}^{NR-1}{d}_{KC}[nr+1,m+1]\mathrm{sinc}\left[\mathsf{\pi }\left(\frac{nr}{KR}-pr\right)\right]{\mathrm{e}}^{-\frac{2\mathsf{\pi }krnr}{KR}\mathrm{j}}}, \\ \forall m\in \overline{0,M-1}, \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}. \end{gathered}$$

(34)

The resulting matrix, $C_{KR,KC}^{m+1}$, can be stored in a 4D array in a cartesian position indicated by the column index $m$ after applying the TDR again. More specifically, if $C$ stands for the 4D array of analysis coefficients, then the matrix $C_{KR,KC}^{m+1}$ is placed in $C$ as follows:

$$C\left[•,•,⌊\frac{m}{KC}⌋+1,m\%KC+1\right]=C_{KR,KC}^{m+1}, \forall m\in \overline{0,M-1},$$

(35)

where “$•$“ stands for “all elements” of the corresponding dimension. Thus, the current element of the 4D block is

$$\begin{gathered} c\left[pr-P{R}_{\min }+1,kr+1,pc-P{C}_{\min }+1,kc+1\right]=c_{KR,KC}^{\left(pc-P{C}_{\min }+1\right)·KC+kc+1}[pr-P{R}_{\min }+1,kr+1], \\ \forall pr\in \overline{P{R}_{\min },P{R}_{\max }}, \forall kr\in \overline{0,KR-1}, \forall pc\in \overline{P{C}_{\min },P{C}_{\max }}, \forall kc\in \overline{0,KC-1}. \end{gathered}$$

(36)

As an alternative, packing the coefficients in a 4D array can be avoided by storing all matrices $C_{KR,KC}^{m+1}$ in a 3D array as layers:

$$C\left[•,•,m+1\right]=C_{KR,KC}^{m+1}, \forall m\in \overline{0,M-1}.$$

(37)

Algorithm 1. Direct JvNT for 2D real-valued signals

3.2. Inverse JvNT Algorithm for Real-Valued 2D Signals

For the synthesis procedure, it is suitable to invoke Algorithm 2 from [3], such that the symmetry property can fully be exploited. Therefore, the previous matrix-oriented approach is useful in the design of this procedure, as well.

Algorithm 2. Inverse JvNT for real-valued 2D signals

4. Simulation Results and Discussion

After implementing the algorithms from Section 3 within the MATLAB™ programming environment, several tests were performed on bw and colored images. All of them are represented by real-valued matrices (bw images) or 3D arrays (colored images). The simulation results for a couple of 2D signals are exhibited and discussed in this section. The JvN dictionaries were configured with an accuracy threshold of $\mathsf{\epsilon }=0.01$.

Two types of tests were performed: a reconstruction test and a redundancy reduction test. The first one targets the image reconstruction accuracy (i.e., the transform invertibility accuracy), depending on the JvNT parameters and the dictionary size. The second one was preferred (among other tests) because there is a direct link between the orthogonality property and the capacity of an orthogonal transform to reduce the redundancy of a signal. Usually, the redundancy is assessed by means of the compression ratio/factor, obtained after selecting the most important transform coefficients, to represent the signal (instead of signals samples). For the 2D JvNT (like in the case of 1D JvNT in [3]), only the theoretical compression capacity was estimated, as the whole compression algorithm depends on various other techniques to be applied on JvNT coefficients afterwards. Thus, the higher the theoretical compression capacity, the stronger the redundancy reduction.

4.1. Black and White Image

The image on top in Figure 2 is a picture of Daisy—a fashionable young singer in the sixties (the photo was downloaded from a public site [31]). In the MATLAB™ programming environment, the photo can be uploaded and displayed as a bw image—see the picture on the left side, on the bottom in Figure 2. Digital bw images like Daisy are represented in gray scale, which means the corresponding matrix includes pixels with integer values varying from 0 (pure black) to 255 (pure white). Each pixel takes one unsigned byte (8 bites) of binary representation. Thus, the visual effect of gray tones is achieved.

Applying JvNT directly on bw image is unwise, as the resulting coefficients are not integers anymore and large errors can be obtained when trying to synthesize the image from them. (In general, for 2D signals with integer values, integer-to-integer transforms are more suitable to employ—see, for example [32,33].) Therefore, the integer pixels were normalized in the range of [−1, +1], with floating point (fp) values. The resulting matrix can be represented as a 3D surface, like that on the right side on the bottom in Figure 2. The image sizes are $NR=658$ and $NC=566$ (with a total of 372,428 pixels).

The 2D FT applied to the normalized image led to the spectrum and phase of Figure 3.

The information of the Daisy image is concentrated in two corners, according to the first window of the figure (on top), where the spectrum on a linear scale is depicted. To the right, the spectrum in dB can be seen. The dominant frequencies are not so obvious under this surface. The phase surface is depicted at the bottom of the figure and reveals almost linear variations along the first frequency axis.

JvN analysis can be performed through Algorithm 1 in two cases:

• for $[KR|KC]=⌊\sqrt{[NR|NC]}+0.5⌋=[26|24]$;
• for $[KR|KC]=⌊\sqrt{2[NR|NC]}+0.5⌋=[36|34]$.

The resulting coefficients cannot be displayed because they belong to 3D or 4D arrays. Note that the total number of coefficients is 4,831,632 for $[KR|KC]=[26,24]$ and 8,029,440 for $[KR|KC]=[36|34]$. Both numbers are very large.

To estimate the theoretical compression capacity of JvNT, the approach is similar to the one in [3]. Thus, the linearized arrays of JvN coefficients were first sorted in descending order of their magnitude. Assume that the total amount of energy accumulated by all coefficients is $\mathsf{E}$ (the sum of squared magnitudes) and set a threshold $\mathsf{\eta }\in [0,1]$ of relative energy. Then, the first coefficients that accumulate an amount of energy at least equal to $\mathsf{\eta }\mathsf{E}$ are counted. Denote their number as $N_{\mathsf{\eta }}\in \mathbb{N}$. In Figure 4, the variations in $N_{\mathsf{\eta }}$ (depending on the threshold $\mathsf{\eta }$) are drawn in both cases. From the figures, one can notice that, if $\mathsf{\eta }=0.95$, then $N_{0.95}=689$ or $N_{0.95}=571$, which represents no more than $[0.014|0.007]\%$ of the total coefficient number, respectively. When comparing to the total number of pixels in the image (i.e., 372,428), both numbers $N_{0.95}$ are very small.

Nevertheless, to be fair, it should be outlined that, based on the difference from the binary representation of pixels, the JvN coefficients are complex-valued and represented in fp double precision, which means that $N_{0.95}$ must be multiplied first by 2, and then, by the number of bytes allocated to double precision, e.g., 8, on a 64-bit computer. In this case, the number of selected bytes is $[\mathrm{11,024}|9136]$, which means approximately $[2.96|2.45]\%$ of the total number of bytes in the image. If the symmetry properties are considered, the number of bytes can be almost halved. Thus, $[1.62|1.35]\%$ of the total number of bytes in the image are sufficient to recover the Daisy picture with fair accuracy.

The two depicted angles standing for the theoretical compression capacity were computed as explained in [3]. The interval $[0,1]$ of $\mathsf{\eta }$ was discretized with a step of 0.01, which resulted in 101 energy thresholds (including the null one). For each couple $\left\{\mathsf{\eta },N_{\mathsf{\eta }}\right\}$, the angle of the first derivative in the origin is estimated as follows:

$${\alpha }_{\mathsf{\eta }}=\arctan \frac{100\mathsf{\eta }}{N_{\mathsf{\eta }}}.$$

(38)

Then, the final angle is obtained by averaging all 100 angles (38) (except for the first one). The higher the average angle, the better the theoretical compression capacity. Figure 4 reveals that the second JvNT, for $[KR|KC]=[36|34]$, is superior to the first one, for $[KR|KC]=[26|24]$, in terms of theoretical compression capacity, since the corresponding average angle is 24.76°, compared to 13.12°.

In

Table 1. Relative energy $\mathsf{\eta }[\%]$ versus JvN coefficient number $N_{\mathsf{\eta }}$ for Daisy image.

$[\mathit{K}\mathit{R}|\mathit{K}\mathit{C}]\downarrow$	$\mathbf{\eta }\to$	20	25	30	35	40	45	50	55	60	65	70	75	80	85	90	95	Angle
$[26|24]$	$N_{\mathsf{\eta }}\to$	73	93	115	137	160	184	208	233	260	288	319	353	389	431	473	689	13.12°
$[36|34]$	$N_{\mathsf{\eta }}\to$	35	45	56	67	79	91	103	116	130	145	162	181	204	233	276	571	24.76°

, the variations in Figure 4 are sampled.

After the completion of the synthesis stage (through Algorithm 2), the obtained results are depicted in Figure 5 (for $[KR|KC]=[26|24]$) and Figure 6 (for $[KR|KC]=[36|34]$).

In both figures, the first synthesized image is obtained in fp, which means the pixels’ normalized values must vary within a [−1, +1] interval. The outliers are saturated, as they are produced by numerical errors. In the upper left corner, one can see the recovered 3D images. Apparently, both images are identical to the original one in Figure 2.

Nevertheless, the recovery errors are non-null, as the windows located in the left corner at bottom of each figure clearly reveals. Although not obvious, the first error (in Figure 5) is slighter smaller, as proven by the performance parameters below. The same behavior is proven by the next column in both figures. This time, the recovered images and errors work with binary (integer) pixels. To obtain the bw image, it suffices to map back the interval [−1, +1] to the integers set at $\overline{0.255}$ by using rounding.

The white and gray spots in the binary error images of figures (see the right-side corners at the bottom) point to the affected pixels (recall that black means a null error). Nevertheless, the recovered bw images (in the upper right corner) are visually identical to the original image on top in Figure 2.

The performance parameters of 2D JvNT are listed in

Table 2. Performance parameters of JvNT in case of Daisy image.

Image Size $\mathit{N}\mathit{R}\times \mathit{N}\mathit{C}$	Sampling Rates $[\mathit{K}\mathit{R}|\mathit{K}\mathit{C}]$	Relative Std of Error [%]			Fitness (Accuracy) [%]			Analysis Runtime [s]	Synthesis Runtime [s]
		fp	Saturated fp	Integer	fp	Saturated fp	Integer
$658\times 566$	$[26|24]$	0.7744	0.643	0.1624	92.81	93.96	97.75	282.74	39.79
	$[36|34]$	0.9574	0.8642	0.2065	91.26	92.05	97.01	455.18	80.84

. The relative standard deviation (std) of the error is computed as follows:

$${\sigma }_{\tilde{A}}=100\frac{{‖A-\tilde{A}‖}_2}{{‖A‖}_2} [\%],$$

(39)

where $A$ is the genuine matrix before analysis, and $\tilde{A}$ is the synthesized matrix. The fitness (which actually evaluates the relative recovery accuracy) is computed by using definitions (19) and (28). More specifically,

$$\mathsf{A}=\frac{100}{1+\frac{{\sigma }_{\tilde{A}}}{10}} [\%].$$

(40)

As one can notice, the recovery accuracy improves while passing from fp representations of pixels to binary representation. The first JvNT is slightly better than the second one in terms of recovering accuracy. Nevertheless, the second JvNT seemingly is better than the first one in terms of compression capacity.

The synthesis runtime is sensibly smaller than the analysis runtime: seven times for $[KR|KC]=[26|24]$ and six times for $[KR|KC]=[36|34]$. The runtime can rapidly increase without using the symmetry properties. Also, in all implementations, the JvN coefficients were recorded in 3D arrays, which also speed up the run (compared to the 4D array solution).

Figure 7 displays the results of lossy synthesis for the Daisy image. In the window on top in the figure, fitness characteristics are drawn for the two JvNTs, depending on the number of strongest coefficients (nsc), $N_{\mathsf{\eta }}$. If the nsc increases, fitness increases, as expected. Although in the beginning (for a small nsc), the second JvNT (for higher sampling rates) performs better in terms of fitness (see the drawing in red), when more and more coefficients are added, the performance of both transforms becomes approximately the same (the first JvNT tacks the lead with very small advancement).

Regarding the fitness characteristics, two points were selected for each couple of sampling rates. Thus, the synthesized images on the left side at the bottom of Figure 7 were obtained for $\mathsf{\eta }=0.95$, i.e., from the strongest $[689|571]$ coefficients (see Table 1 again). Both images are visibly distorted (blurred) because the number of selected coefficients is too small (despite them storing approximately 95% of image energy). This is confirmed either by the recovery errors with a relative std of $[4.61|4.51]\%$ or by the fitness values of $[68.45|68.93]\%$. When closely looking at the images, one can see that the second JvNT (at the bottom) offered slightly better accuracy (and slightly smaller error to the right) than the first JvNT (on top) (in fact, the resulting fitness is about $0.5\%$ higher).

If we move now to the right on fitness characteristics, up to the points determined by $\mathsf{\eta }=0.99$, this time, the nscs are $[4328|4815]$. In the 64-bit representation, they take $[\mathrm{138,496}|\mathrm{154,080}]$ bytes, i.e., $[37.19|41.37]\%$ of the image’s total number of bytes (372,428). With the symmetry property, these percentages are almost halved: $[19.5|21.8]\%$ of the image’s total number of bytes. This yields sensibly, reducing the blurring effect in both synthesized images on the right side at the bottom of the figure. The recovery errors are smaller, as well, with a relative std of $[1.61|1.59]\%.$ Consequently, fitness increases to $[86.1|86.27]\%$. In this case, too, the second JvNT performs slightly better in terms of fitness, but at the expense of a higher nsc employed in the synthesis.

4.2. Colored Image

The data type required to operate with colored images using a computer is 3D array. If the image is represented in a red–green–blue (RGB) system (which is one of the most employed), then the array consists of three layers—one for each color component. On every layer, unsigned integers determine the proportion of each color that contributes to pixel definition, as suggested in Figure 8.

So, practically, a colored image is represented by three matrices of bytes. The image under study in this article is depicted in Figure 9: a photo of the masterpiece painting Ship of Dreams (also known as The Wind) by Salvador Dali (the photo was taken from a private paintings album).

The image sizes are $NR=540$ and $NC=720$ (with a total of 388,800 pixels represented on three unsigned bytes each, i.e., 1,166,400 bytes for the whole 3D array). Apparently, beside the colormap, the layers are not so different from each other, as all three must follow the same pattern imposed by the painting itself. Nevertheless, the three matrices are different, as proven by their norms: 58,883.81 for the R layer; 83,152.92 for the G layer; 89,721.65 for the B layer. Also, the norms computed on the differences between them are quite big: 32,776.94 for R-G; 45,320.43 for R-B; 16,823.66 for G-B. Interestingly, the G and B layers are closer to each other than to the R layer. Moreover, their norms are bigger. These numerical results are visually confirmed by the whole image, where one can easily notice that the blue and green colors are dominant over the red color, which only fills small areas.

Two approaches can be considered when applying JvNT to colored images: extend definition (12) to 3D arrays (i.e., work with three instead of two dictionaries) or process every layer separately by means of 2D transform. The first approach is appropriate when the number of layers is large enough. In the case of colored images, with only three layers, the third dictionary would be devoted to processing 1D signals with a length of three, which is unnecessarily complex. Therefore, the second approach (process every layer as a 2D signal) is more suitable. This strategy is like the one employed for the bw Daisy image (first, normalize the unsigned integers of the layer; then, process the resulting matrix with the fp values; finally, come back to the unsigned integers for the recovered layer).

The 2D FT applied on each layer led to the results in Figure 10.

The spectral representation of all layers has the same main characteristics as for the bw image: the energy/information seems to be concentrated in the main two corners of the frequency plane, although there is energy dissipation on the remaining Fourier coefficients. Therefore, compression of the image by using FT could perform as well.

Algorithm 1 was run for each layer in the case of the sampling rates below:

• $[KR|KC]=⌊\sqrt{[NR|NC]}+0.5⌋=[23|27]$;
• $[KR|KC]=⌊\sqrt{2[NR|NC]}+0.5⌋=[33|38]$.

The total number of analysis coefficients is quite big: 4,862,430 per layer for $[KR|KC]=[23|27]$ and 8,226,240 per layer for $[KR|KC]=[33|38]$. Nevertheless, like in case of Daisy photo, only a few of them encode essential information about the analyzed image. A similar analysis of compression capacity led to the results in Figure 11 and Figure 12 and

Table 3. Relative energy $\mathsf{\eta }[\%]$ versus JvN coefficient number $N_{\mathsf{\eta }}$ for Ship of Dreams/The Wind painting layers.

$[\mathit{K}\mathit{R}|\mathit{K}\mathit{C}]\downarrow$	$\mathbf{\eta }\to$	20	25	30	35	40	45	50	55	60	65	70	75	80	85	90	95	Angle	Layer
$[23|27]$	$N_{\mathsf{\eta }}\to$	54	70	87	105	125	146	169	194	221	253	292	340	410	574	1199	4082	15.28°	R
		36	47	61	75	92	111	133	160	193	236	294	390	605	1139	2623	8707	19.14°	G
		55	73	91	112	133	157	182	209	240	274	315	367	436	543	1031	4059	14.64°	B
$[33|38]$	$N_{\mathsf{\eta }}\to$	27	35	44	54	64	76	88	102	118	137	162	204	297	526	1185	4496	26.88°	R
		19	26	33	42	52	63	77	93	114	143	192	295	526	1083	2716	9328	29.80°	G
		29	38	48	58	70	82	95	109	125	145	170	203	250	362	887	4113	25.49°	B

.

If $\mathsf{\eta }=0.95$, the total number of coefficients (on all three layers) is 16,848 for $[KR|KC]=[23|27]$ and 17,937 for $[KR|KC]=[33|38]$. This means that 95% of the image energy is concentrated in $[0.1155|0.0727]\%$ of the total coefficient number. In the case of 64-bit representation, the number of bytes that can be selected is 269,568 for $[KR|KC]=[23|27]$ and 286,992 for $[KR|KC]=[33|38]$, i.e., $[23.11|24.61]\%$ of the total number of bytes in the image (on all three layers).

If the symmetry properties are considered, the percentages are almost twice smaller: $[12.66|13.51]\%$. The average on the three layers of the compression angles depicted in Figure 11 is 16.35° (for $[KR|KC]=[23|27]$), while this average increases to 27.39° for the angles in Figure 12 (where $[KR|KC]=[33|38]$). This result shows that, similarly to the analyzed bw image, the second JvNT (for higher sampling rates) seemingly can better compress the color image than the JvNT for smaller sampling rates.

Each layer was recovered from the corresponding JvN coefficients by using Algorithm 2. The synthesis results are displayed in Figure 13 and Figure 14 for $[KR|KC]=[23|27]$ and $[KR|KC]=[33|38]$, respectively. In both figures, the layers are shown in rows (red on top, green in the middle and blue at the bottom). The synthesized paintings are recomposed from RGB layers like in Figure 15.

When looking at all the previous pictures, one can straightforwardly conclude that the synthesized images are very close to the original ones. This subjective observation is confirmed by the performance parameters in

Table 4. Performance parameters of JvNT in case of Ship of Dreams/The Wind painting.

Image Size $\mathit{N}\mathit{R}\times \mathit{N}\mathit{C}$	Sampling Rates $[\mathit{K}\mathit{R}|\mathit{K}\mathit{C}]$	Layer	Relative Std of Error [%]			Accuracy [%]			Analysis Runtime [s]	Synthesis Runtime [s]
			fp	Saturated fp	Integer	fp	Saturated fp	Integer
$540\times 720$	$[23|27]$	R	0.8122	0.8032	0.5187	92.49	92.57	95.07	263.05	36.14
		G	1.1138	1.1121	0.3266	89.98	89.99	96.84	268.28	37.24
		B	0.9174	0.9024	0.3335	91.60	91.72	96.77	269.43	36.21
		P	0.9478	0.9393	0.3929	91.36	91.43	96.23	800.67	109.59
	$[33|38]$	R	0.9030	0.8953	0.5795	91.72	91.78	94.52	425.91	73.15
		G	1.0567	1.0552	0.3103	90.44	90.46	96.99	436.25	80.64
		B	0.8057	0.7900	0.2907	92.54	92.68	97.18	433.21	80.96
		P	0.9218	0.9135	0.3935	91.57	91.64	96.23	1295.37	234.75

.

The notation “P” in the table above stands for the whole picture, as composed by all three layers. On each row “P”, the averages of the errors relative the std and fitness values are listed. Also, the runtimes are cumulated on row “P”, as if the layers were sequentially processed. One can see that the two JvNTs are similar in terms of recovery accuracy, although the second one runs longer (however, recall that this JvNT seems to perform better in terms of compression capacity). The synthesis runtime is 6 to 8 times smaller than the analysis runtime, especially thanks to the symmetry properties.

The analysis of lossy synthesis was performed like in the case of the Daisy image. Nevertheless, for Dali’s painting, each layer was processed separately. The obtained results are illustrated in Figure 16 (for $[KR|KC]=[23|27]$) and Figure 17 (for $[KR|KC]=[33|38]$).

In both figures, the synthesis results are depicted at the bottom, for $\mathsf{\eta }=0.95$ to the left and for $\mathsf{\eta }=0.99$ to the right. In both cases, the blurring effect is obvious for the smaller threshold. In turn, when the threshold increases, the synthesized images are clearer. Below the images, the recovery errors are shown for each layer (red, green and blue from left to right) after passing to integer pixels. The same representation scale was employed for both thresholds. One can see that the magnitude of errors for $\mathsf{\eta }=0.99$ is smaller, as expected. The fitness variation is drawn above recovered images for each layer. A fourth variation is added, standing for the whole picture, by averaging the layer variations. Unlike for the bw image, in this case, the threshold $\mathsf{\eta }$ led to different nscs, depending on each layer. In

Table 5. Number of strongest coefficients and fitness (accuracy) values in case of Ship of Dreams/The Wind painting.

Sampling Rates $[\mathit{K}\mathit{R}|\mathit{K}\mathit{C}]$	Layer	$\mathbf{\eta }=0.95$		$\mathbf{\eta }=0.99$
		${\mathit{N}}_{0.95}$	Accuracy [%]	${\mathit{N}}_{0.99}$	Accuracy [%]
$[23|27]$	R	4082	73.48	30,814	85.67
	G	8707	84.43	47,479	91.16
	B	4059	80.98	32,046	90.86
	P	16,848	79.63	110,339	89.23
$[33|38]$	R	4496	71.58	35,823	83.98
	G	9328	84.57	53,582	91.56
	B	4113	81.24	34,458	91.18
	P	17,937	79.13	123,863	88.91

, the number of such coefficients, together with the corresponding fitness values, are listed for the two outlined groups of points in each figure.

The Green layer seems to be the most accurately recovered, but at the expense of a greater number of employed coefficients (almost twice bigger than for the other layers). This means the Green layer is less auto-correlated than the other two layers. For the whole picture, the nscs were cumulated, while the accuracy was averaged. Although the nsc is quite small for $\mathsf{\eta }=0.95$ (as analyzed above), the synthesized images have modest visual quality. In the case of $\mathsf{\eta }=0.99$, the visual quality is improved. This time, though, the nscs are 6–7 times bigger. More specifically, the nscs are $[\mathrm{110,339}|\mathrm{123,863}]$. They take $[0.7564|0.5019]\%$ of the total coefficient number. In 64-bit representation, the numbers of bytes that correspond to $[\mathrm{110,339}|\mathrm{123,863}]$ are $[\mathrm{1,765,424}|\mathrm{1,981,808}]$, respectively (recall that JvN coefficients are complex-valued). Obviously, both numbers are bigger than the number of bytes of the whole image (1,166,400). Employing the symmetry property becomes crucial for reducing the total number of bytes. In this case, the image can be recovered with good visual quality from approximately $[80.68|89.95]\%$ of the total number of bytes in the image (on all three layers).

The explanation for the last result compared to the previous one, coming from the Daisy image (where the theoretical compression rate was much better), very likely resides in the fact that Dali’s painting reveals a much more complex fractal structure. JvNT is a smooth transform with limitations in analyzing and/or the compression of highly fractal images. Normally, the analysis instrument should have the same nature as the analyzed object. In this case, fractal transform (for example, built by using Daubechies’ wavelets [34]) could perform better. Dali’s painting was selected on purpose to outline, on one hand, that symmetry becomes crucial for transform effectiveness and, on the other hand, that, like any other transform, JvNT has limitations, too, in this case, imposed by its non-fractal nature.

The analysis of lossy synthesis is just the starting point for achieving a sounder analysis concerning the lossy compression of images. When computing the compression ratio, the side information must be accounted for, as well. This means either encoding and sending the positions of the strongest coefficients or sending all coefficients, after setting to null (or even to smaller values) the weakest coefficients. In both cases, the compression performance is reduced, as the total number of bytes necessary to recover the image increases. Usually, after some orthogonal transform is applied to the image, a series of compression algorithms are employed to encode the resulting coefficients. This is a sophisticated endeavor, beyond the goal of the current article.

5. Concluding Remarks

The extension of John von Neumann’s orthogonal transform from time-varying (1D) signals to space-varying (2D) signals was the main goal of this article. The requirement to meet was to preserve both the orthogonality and the invertibility properties of resulted 2D transform. Two approaches were investigated in this study: a rotational one and a cartesian one. Because, in the first approach, the orthogonality is very difficult to preserve, the second approach was adopted. Thus, the natural link with 1D transform was exploited, on one hand, to keep intact the orthogonality property and, on the other hand, to design efficient analysis–synthesis algorithms for image processing. Moreover, it has been shown that using a couple of 1D transforms for 2D signals is equivalent to the analysis and synthesis of such signals with 2D transform, thanks to a fortunate factorization of the cardinal sinus. After running the corresponding numerical algorithms on black and white as well as on colored images, the results were analyzed in two main respects: the theoretical compression capacity and the quality of the reconstructed image. According to the obtained results, one can conclude that John von Neumann’s 2D transform is a promising instrument to be applied to images prior to other compression methods. Beside compression, 2D transform can be employed in the analysis of any 2D signal (not necessarily an image) and, furthermore, the extension of its definition to multi-dimension signals is straightforward, as, for each new dimension, one new 1D transform can be defined and employed.