Introduction

Hyperspectral imaging which combines spectroscopy and digital imaging has become one of the powerful tools for solving many computer vision tasks such as quality inspection and sorting efficiently. Hyperspectral imaging has been widely used in remote sensing applications by virtue of modern sensor technologies which can cover large surfaces of earth. In recent years, ground-based hyperspectral imaging has emerged as a powerful tool for quality and safety evaluation of food items, forensic science, medical surgery and diagnosis, military applications, and restoration of artworks [1]. HSI is more effective in detecting the external or internal quality and chemical compositions of horticultural products [2] and therefore there is a tremendous increase in the research on hyperspectral imaging used for quality and safety checking of fruits and vegetables. A hyperspectral imaging system obtains a two dimensional (2D) matrix for each wavelength ranging from VIS (visible) to NIR (Near Infrared). The resulting structure is a three dimensional (3D) image dataset, which is called a hypercube. The data acquisition of hyperspectral images can be done in 4 modes-line scanning, area scanning, point scanning, and single shot [2]. In a pushbroom based hyperspectral system, the sensor captures spectral information at each line, and by moving the camera or the object, the system captures the entire spatial region [2]. The major components of a pushbroom hyperspectral imaging system are: illumination, camera, objective lens, transportation plate, and a data processing unit, as shown in Fig. 1 [3]. The spectral range produced by the system depends on illumination conditions, camera, and the typical range of a VNIR system is 400 to 1000 nm. Even if HSI is a significant The anomalous observations can be due to several reasons, such as the dysfunction and noise in the sensor or from different stages in the workflow, shape and geometry of the scene, and the radiation technique [4]. In this paper, we try to develop a new model for removing specific noise artefacts affecting the HS images.

The pushbroom sensor used for image acquisition produces stripes, and random noise [5]. The random noise can be viewed as an adaptive white Gaussian(AWGN) random process with variance 𝜎 2 and mean zero. Dead detectors will give rise to single line drop-out for pushbroom sensors [6]. The neighbouring elements in the CCD(Charge Coupled Device) array may have sensitivity variations, which causes vertical stripe noise [7]. HSIs also suffer from impulse noise [8]. These noises can affect the subsequent internal and external quality evaluation and defect detection processes. Improving the HSI quality merely through a hardware scheme is unsustainable and impractical. Therefore, suitable image pre-processing techniques are to be applied to obtain a high-quality HSI before subsequent applications.

HS imaging was actually developed for remote sensing applications, and so far, a huge number of proposals have been made in the HSI denoising. These techniques can be broadly classified into five categories [9]. In the first category, each band in the HSI is treated as an image, and then commonly used two dimensional (2D) denoising methods are applied in each band to remove noise [10], [11]. But the results are not promising since the band-wise processing does not consider the strong correlations between bands. Another method is to treat the HSI as a multidimensional data cube, and the volumetric data denoising methods can accomplish the denoising [12], [13]. However, the results are unsatisfactory because these methods do not consider the correlation among the bands. Therefore, it is necessary to consider the spatial and spectral information simultaneously to enhance the denoising quality. The third type of HSI denoising method combines spatial and spectral information, and several algorithms have been proposed based on this concept.e.g. [14], [15]. The results obtained from this fusion method are promising compared with single view(either spatial or spectral) methods. A few transform domain-based techniques were also proposed in this category [16], [17].

The fourth type of HSI denoising method falls under the category of low rank based techniques. In [18], a randomized singular value decomposition is introduced and a noise adjusted iterative framework is proposed for low rank matrix approximation (NAILRMA). A SSTV based low Tucker rank model (LRTDTV) was proposed to preserve the spectral signatures of HSI in [19]. The method proposed by Chen et al. [20] is based on weighted group sparsity regularization, and this concept is incorporated for low rank tucker decomposition (LRTDGS), in which results were shown to be improved in comparison with the previous TV based methods. Zheng et al. [21] proposed the fibered rank and suggested the convex (3DTNN) and non-convex models (3DLogTNN) to separate clean HSI from noisy images.

Recently, many deep learning based algorithms have been proposed, and their performance is encouraging [22], [23]. But their results are highly dependent on the quality and volume of training data which limits their application.

Restoration of HS images for remote sensing applications has been well studied in the literature. But the researchers have less focused on restoring the images for close range applications. Moreover, the algorithms used for remote sensing applications may not be suitable for HSI restoration of data captured at proximity range. This is evident from Fig. 2 The low rank based HSI restoration techniques can be broadly classified as matrix based and tensor based approaches. Even though several methods have been proposed based on low rank matrix modelling for removing mixed noise, these methods failed to utilise the correlation in spatial and spectral modes effectively, leading to suboptimal denoising results under severe noisy conditions. So many recent research works included the direct tensor modelling techniques, and it has been found that these techniques are superior to matrix based techniques in computing higher order data. Such studies motivated us to propose a new optimisation model for restoring noise-contaminated hyperspectral images using direct tensor modelling techniques exploiting the spatial and spectral information.

Our Contributions:

. As per our knowledge, this is the first effort to realise a tensor based low rank(LR) and SSTV model for the removal of noise artefacts from the HSI data captured at proximity range.

. Further, 𝑙 1 norm and Frobenius norms are used to address the sparse noise and heavy Gaussian noise in the images.

. The formulated optimisation problem is decomposed into several subproblems and solved using the Alternating Direction Method of Multipliers (ADMM).

. The algorithm is implemented successfully on the fruit HSI data with simulated noise.

The rest of the paper is organised as follows. Section 2 explains the tensor based low rank and SSTV model and the optimisation solution in detail. The experimental results are presented in section 3. Finally, the conclusions are derived in section 4.

Proposed Model

A noisy HSI cube ℋ may be represented as,

ℋ = 𝒥 + ℰ(1)

where 𝒥 represents the clean HS image and ℰ represents noise artefacts. The aim of HSI restoration task is to obtain the noise-free HS image 𝒥 from the noise-contaminated image ℋ.

The noise term includes Gaussian noise 𝒲 and the sparse noise 𝒮. Dead pixels, stripes and impulse noise fall into the category of sparse noise. Based on this, the degraded HSI can be modelled as,

ℋ = 𝒥 + 𝒮 + 𝒲(2)

Using the TRPCA (Tensor Robust Principal Component Analysis) model [24], the image can be expressed as a combination of a rank function which describes the low rank property and a ℓ 0 norm to represent the sparse noise.

argmin 𝒥 rank(𝒥 ) + 𝜆‖𝒮‖ 0 𝑠.𝑡. ℋ = 𝒥 + 𝒮(3)

The optimization problem Eq.( 3) is a nonconvex optimization problem. So we use the convex surrogates, tensor nuclear norm ‖.‖ ⊛ and ℓ 1 norm ‖.‖ 1 to approximate the rank(𝒥 ) and ‖𝒮‖ 0 , respectively. The reformulated convex optimization problem is, argmin

𝒥 ‖𝒥 ‖ ⊛ + 𝜆‖𝒮‖ 1 + 𝛽‖𝒲‖ 2 F , 𝑠.𝑡. ℋ = 𝒥 + 𝒮 + 𝒲(4)

where ℋ,𝒥 ,𝒮, and 𝒲 are 3 𝑟𝑑 order tensors. The low rank approximation can be obtained using tensor singular value decomposition algorithm [25]. The objective of using the TV regularization is to take advantage of the piecewise smoothness property in both the spectral and spatial domains. The widely used 2D TV regularizer consider only the spatial smoothness structure of HSI. So we also use a SSTV regularizer to explore the spectral smoothness.

‖𝒥 ‖ SSTV := ∑︁ i,j,k w 1 |p i,j,k − p i,j,k-1 | + w 2 |p i,j,k − p i,j-1,k | + w 3 |p i,j,k − p i-1,j,k |(5)

where 𝑝 𝑖,𝑗,𝑘 , is the (i, j, k) 𝑡ℎ entry of 𝒥 . The term 𝑤 𝑙 (l = 1, 2, 3) is the weight along the l 𝑡ℎ mode of 𝒥 . These weights decide the strength of regularization. When the SSTV regularization is also added to the rank-constrained RPCA, the model will become, argmin

𝒥 ‖𝒥 ‖ ⊛ + 𝜆‖𝒮‖ 1 + 𝜏 ‖𝒥 ‖ SSTV , 𝑠.𝑡. ℋ = 𝒥 + 𝒮 + 𝒲(6)

By adding some auxiliary variables Eq.( 6) can be rewritten as,

argmin 𝒥 ‖𝒥 ‖ ⊛ + 𝜆‖𝒮‖ 1 + 𝛽‖𝒲‖ 2 F + 𝜏 ‖ℱ‖ 1 , 𝑠.𝑡. ℋ = 𝒥 + 𝒮 + 𝒲, 𝒥 = 𝒵, 𝐷 𝑤 (𝒵) = ℱ(7)

where 𝐷 𝑤 (•) is the third order weighted difference operator. It can be given by

𝐷 𝑤 (•) = [ 𝑤 𝑚 × 𝐷 𝑚 (•); 𝑤 𝑛 × 𝐷 𝑛 (•); 𝑤 𝑘 × 𝐷 𝑘 (•)]

where 𝐷 𝑚 , 𝐷 𝑛 , and 𝐷 𝑘 represents the first order difference operators in the different orientations of a 3D HSI. The optimization problem in Eq.( 7) can be rewritten using the augmented Lagrangian multiplier (ALM) method:

argmin 𝒥 ‖𝒥 ‖ ⊛ + 𝜆‖𝒮‖ 1 + 𝛽‖𝒲‖ 2 F + 𝜏 ‖ℱ‖ 1 + ⟨𝐺 1 , ℋ − 𝒥 − 𝒮 − 𝒲⟩ + 𝜇 2 ‖ℋ − 𝒥 − 𝒮 − 𝒲‖ 2 F + ⟨𝐺 2 , 𝒥 − 𝒵⟩ + 𝜇 2 ‖𝒥 − 𝒵‖ 2 F + ⟨𝐺 3 , 𝐷 𝑤 (𝒵) − ℱ ⟩ + 𝜇 2 ‖𝐷 𝑤 (𝒵) − ℱ ‖ 2 F (8)

𝒥 update

𝒥 𝑘+1 =argmin 𝒥 ‖𝒥 ‖ ⊛ + ⟨𝐺 k 1 , ℋ − 𝒥 − 𝒮 k − 𝒲 k ⟩ + 𝜇 2 ‖ℋ − 𝒥 − 𝒮 k − 𝒲 k ‖ 2 F + ⟨𝐺 2 , 𝒥 − 𝒵 k ⟩ + 𝜇 2 ‖𝒥 − 𝒵 k ‖ 2 F =argmin 𝒥 ‖𝒥 ‖ ⊛ + 2( 𝜇 2 )‖𝒥 − ℳ‖ 2 F (9)

where

ℳ = 1 2 (︂ ℋ + 𝒵 k − 𝒮 k − 𝒲 k + (𝐺 k 1 − 𝐺 k 2 ) 𝜇 )︂ 𝒮 update 𝒮 𝑘+1 =argmin 𝒮 𝜆‖𝒮‖ 1 + ⟨𝐺 k 1 , ℋ − 𝒥 k+1 − 𝒮 − 𝒲 k ⟩ + 𝜇 2 ‖ℋ − 𝒥 k+1 − 𝒮 − 𝒲 k ‖ 2 F =argmin 𝒮 𝜆‖𝒮‖ 1 + 𝜇 2 ‖𝒮 − (ℋ − 𝒥 k+1 − 𝒮 − 𝒲 k + 𝐺 1 𝜇 )‖ 2 F =𝑠𝑜𝑓 𝑡 𝜆 𝜇 (︂ ℋ − 𝒥 k+1 − 𝒲 k + 𝐺 k 1 𝜇 )︂(10)

𝒲 update

𝒲 𝑘+1 =argmin 𝒲 𝛽‖𝒲‖ 2 F + ⟨𝐺 𝑘 1 , ℋ − 𝒥 k+1 − 𝒮 k+1 − 𝒲⟩ + 𝜇 2 ‖ℋ − 𝒥 k+1 − 𝒮 k+1 − 𝒲‖ 2 F =argmin 𝒲 (𝛽 + 𝜇 2 ) ⃦ ⃦ ⃦ ⃦ 𝒲 − 𝜇(ℋ − 𝒥 k+1 − 𝒮 k+1 ) + 𝐺 k 1 𝜇 + 2𝛽 ⃦ ⃦ ⃦ ⃦ 2 F = 𝜇(ℋ − 𝒥 k+1 − 𝒮 k+1 ) + 𝐺 k 1 𝜇 + 2𝛽(11)

𝒵 update

𝒵 𝑘+1 =argmin 𝒵 ⟨𝐺 k 2 , 𝒥 k+1 − 𝒵⟩ + 𝜇 2 ‖𝒥 k+1 − 𝒵‖ 2 F + ⟨𝐺 k 3 , 𝐷 𝑤 (𝒵) − ℱ k ⟩ + 𝜇 2 ‖𝐷 𝑤 (𝒵) − 𝐹 k ‖ 2 F (12)

The solution to the above equation can be obtained from the following linear equation.

(𝜇I + 𝜇𝐷 𝑇 𝑤 𝐷 𝑤 )𝒵 = 𝜇𝒥 𝑘+1 + 𝜇𝐷 𝑇 𝑤 (ℱ 𝑘 ) + 𝐺 𝑘 2 − 𝐷 𝑇 𝑤 (𝐺 𝑘 3 )

where 𝐷 𝑇 𝑤 indicates the adjoint operator of 𝐷 𝑤 . We adopt the Fourier Transform to solve Eq. (12).

H z = 𝜇𝒥 k+1 + 𝜇𝐷 𝑇 𝑤 (ℱ 𝑘 ) + 𝐺 𝑘 2 − 𝐷 𝑇 𝑤 (𝐺 𝑘 3 ) Tz = 𝑤 2 𝑚 |fftn(𝐷 𝑚 )| 2 + 𝑤 2 𝑛 |fftn(𝐷 𝑛 )| 2 + 𝑤 2 𝑘 |fftn(𝐷 𝑘 )| 2 𝒵 k+1 = ifftn (︂ fftn(𝐻 𝑧 ) 𝜇I + 𝜇𝑇 𝑧 )︂(13)

ℱ update

ℱ 𝑘+1 =argmin ℱ 𝜏 ‖ℱ‖ 1 + ⟨𝐺 k 3 , 𝐷 𝑤 (𝒵 k+1 ) − ℱ ⟩ + 𝜇 2 ‖𝐷 𝑤 (𝒵 k+1 ) − ℱ ‖ 2 F =argmin ℱ 𝜏 ‖ℱ‖ 1 + 𝜇 2 ‖ℱ − (𝐷 𝑤 (𝒵 k+1 ) + 𝐺 k 3 𝜇 )‖ 2 F =𝑠𝑜𝑓 𝑡 𝜏 𝜇 (𝐷 𝑤 (𝒵 𝑘+1 ) + 𝐺 𝑘 3 𝜇 )(14)

The updates of the multipliers are given by 9),( 10),( 11),( 12),( 13) and ( 14)

𝐺 k+1 1 = 𝐺 k 1 + 𝜇(ℋ − 𝒥 k+1 − 𝒮 k+1 − 𝒲 k+1 ) 𝐺 k+1 2 = 𝐺 k 2 + 𝜇(𝒥 k+1 − 𝒵 k+1 ) 𝐺 k+1 3 = 𝐺 k 3 + 𝜇(𝐷(𝒵 k+1 ) − ℱ k+1 )(157

update 𝐺1, 𝐺2, 𝐺3 via Eq. ( 15)

8

update parameter 𝜇 :=min(𝜌𝜇, 𝜇𝑚𝑎𝑥)

9

check for the convergence criteria

10 ‖𝒥 k −𝒥 k+1 ‖ 2 𝐹 ‖ℋ‖ 2 𝐹 ≤ 𝜖 11 end

EXPERIMENTAL RESULTS AND DISCUSSION

The proposed method is compared with five other best-performing denoising techniques to show the effectiveness in HSI restoration. These methods include LLRGTV [26], 3DTNN [21], 3DLogTNN [21], NAILRMA [18], and LRTDTV [19]. The experiments for all the compared methods are conducted using the same parameters given by the authors in the respective papers to ensure optimum performance. An Intel Xeon CPU@ 3.50 GHz with 64-GB RAM is used to conduct the experiments in MATLAB R2022a

Dataset and Experiments:

Hyperspectral images with different spectral ranges can be used to detect the quality attributes and defects of horticultural products [27]. The dataset for this study is obtained using a camera HySpex-VNIR-1800 which has a spectral sensitivity from 400 nm to 1000 nm and spectral sampling of 3.18 nm [3]. This VNIR push broom scanner records 186 spectral bands and 1800 pixels across the field of view. It uses a polariser to avoid specular reflection. In the present case study, we used the fruits cherry and strawberry. Since the size of the dataset is computationally intensive, it is cropped to 500×500 pixels for cherry and 600×700 pixels for strawberry. Real HSIs are affected by a variety of noises, such as Gaussian noise and sparse noise in different amounts. In order to simulate real noise conditions, we consider different combinations of these noises in varying proportions. Thus we can define the following noise cases as follows.

Noise case 1: Gaussian noise + Impulse noise: Gaussian noise and salt and pepper noise with equal distribution are added to each and every band. The variance of the Gaussian noise is set as 0.2 with zero mean, and the percentage of the impulse noise added is 0.2.

Noise case 2: Gaussian noise + Impulse noise + Dead lines: Here, the variance of Gaussian noise and percentage of impulse noise were fixed as 0.15. Besides this, we add dead lines whose width varies randomly between 1 and 3. The dead lines are added to bands ranging from 101 to 140, and the number of lines randomly varies between 3 and 10.

Noise case 3: Gaussian noise + Impulse noise + Dead lines + Stripes: Here, in addition to the noises as in case 2, we add some stripes also to bands ranging from 91 to 130. The number of stripes randomly varies from 30 to 40. The variance of the Gaussian noise was 0.05 with zero mean, and the percentage of the impulse noise was 0.05.

The visual quality and quantitative metrics for the above noise cases are tested for all methods used for comparison. Visual quality comparison: The denoising results for case 1 and case 3, of band 120 are shown in Figs. 3,4,5 and 6. From Figs. 3(b),4(b),5(b) and 6(b) it can be observed that the original HSIs are severely affected by a combination of various noises. In all the noise cases, it is found that the proposed method gives the best performance. This demonstrates that the low rank minimisation, together with the TV regularisation, perform well for HSI restoration. Quantitative comparison: For quantitative evaluation, we use the metrics such as the mean peak signal-to-noise ratio (MPSNR) and the mean structural similarity index (MSSIM). The higher the values for the PSNR and SSIM, the better the restoration results.

MPSNR = 1 𝑝 𝑝 ∑︁ 𝑗=0 𝑝𝑠𝑛𝑟 𝑗 MSSIM = 1 𝑝 𝑝 ∑︁ 𝑗=0 𝑠𝑠𝑖𝑚 𝑗

where 𝑝𝑠𝑛𝑟 𝑗 and 𝑠𝑠𝑖𝑚 𝑗 are the PSNR and SSIM values for the 𝑗 th band, respectively. Table 1 and Table 2 show the restoration results of cherry and strawberry respectively for all methods considered for comparison. The psnr and ssim values obtained for the different noise cases clearly shows that the proposed method gives better restoration results than the other methods.

Discussion:

Parameter setting: There are several parameters in the algorithm that need to be carefully adjusted in order to get optimum restoration results. The sparsity regularization parameter 𝜆 can be selected as 𝜆 = 100x 𝐶 √ 𝑎𝑏 , where a is the height and b is the width of a single HSI band, and C is a parameter for tuning. In all the simulated data experiments, the value of C is fixed as 10. 𝛽 is the Frobenius norm regularization parameter, which is selected as the inverse of the variance of Gaussian noise. The SSTV regularization parameter 𝜏 can be fixed as a constant equal to 1. The weights of the SSTV are chosen as 1 along the spatial mode and a value in the range of 0 to 1 in the spectral mode. The typical value used in the simulation experiment is data . Table 3 and Table 4 show the details for different methods. It can be observed that the NAILRMA method has the lowest computation time among all the compared methods. But, the quantitative and visual comparison results of this method are not good. The converging time for LLRGTV, 3DTNN and 3DLogTNN methods are much higher than the other methods since they require running an expensive iterative algorithm. The computation time of the proposed method is relatively higher than some other methods such as LRTDTV and NAILRMA but significantly lower than LLRGTV, 3DTNN and 3DLogTNN.

Conclusion

A low rank tensor based model is proposed for removing noise artefacts in the line scan-based hyperspectral images. Additionally an SSTV regularization is used, which preserves the spatial smoothness and spectral correlation. We have used 𝑙 1 norm to detect the sparse noise effectively, and in order to tackle the heavy Gaussian noise conditions which may occur in real situations, we have also considered a Frobenius norm. A new algorithm based on the ALM method is designed to solve the resulting nonconvex optimization model. The algorithm is implemented and tested successfully on the fruit HSI data with simulated noise. The experiments show that the proposed method performs better than some popular methods used for comparison in terms