Exploring the potential of NIR hyperspectral imaging for automated 1 quantification of rind amount in grated Parmigiano Reggiano cheese

Manuscript


Introduction
Parmigiano Reggiano (P-R) is a long-ripened, cooked, hard cheese produced in Italy and registered with Protected Denomination of Origin (PDO).P-R represents one of the most important typical Italian food products and it is exported worldwide.This cheese is manufactured from raw and unheated bovine milk, and the whole production chain must take place in a restricted area in Northern Italy (Malacarne et al., 2008).
The PDO is extended also to grated cheese obtained from Parmigiano Reggiano cheese wheels, provided that the product is grated in the specific production area and packaged immediately afterwards, in order to avoid modifications of its organoleptic properties.
Furthermore, grated cheese products designated as Parmigiano Reggiano should meet technical and technological parameters ruled by the Specifications of Parmigiano Reggiano Cheese (https://www.parmigianoreggiano.com/consortium/rules_regulation_2/default.aspx ), that regulates all stages of P-R production, including cow feeding, cheese manufacturing and ripening process.
One of the different quality parameters of grated P-R cheese regulated by the Specifications is the amount of rind.The rind is the external part of cheese wheels; although edible, in long-ripened cheeses it has chemical and physical properties different from those of the inner part of the cheese.
These differences are mainly caused by exposure to environmental conditions during ripening, which determines a decrease in moisture content, proteolytic activity and a higher degree of oxidation (Cattaneo et al., 2008;Karoui et al., 2007;Malacarne et al., 2019).Furthermore, the different properties of rind and pulp also affect size and shape of grated particles.In fact, generally rind particles are finer and less circular than those derived from the pulp (Alinovi et al., 2019).As a consequence of its peculiar chemical and textural properties, an excessive amount of rind in grated cheese can be perceived by consumers, negatively affecting the organoleptic properties of the product (Zannoni and Hunter, 2015).For these reasons, the percentage of rind in grated P-R cheese should not exceed the 18% (w/w) threshold value.
In order to ensure quality compliance of commercial products of grated P-R cheese and to avoid counterfeits, it is essential to implement effective analytical methods to correctly estimate the rind amount, meeting the requirements of low costs, limited sample preparation and short times of analysis.In this context, Near Infrared (NIR) spectroscopy has been widely employed for fast and non-destructive analysis and characterization of food products, thanks to its ability to easily provide a spectral fingerprint codifying the chemical composition of the analysed sample (Curda et al., 2004;Woodcok et al., 2008;Foca et al., 2013;Kraggerud et al.;2014).In particular, previous research studies evaluated the effectiveness of NIR spectroscopy in verifying the authenticity of P-R grated cheese and in discriminating compliant from non-compliant samples (Cevoli et al., 2013a, Cevoli et al., 2013b).However, grated cheese is an inhomogeneous food matrix composed by unevenly dispersed particles derived from both cheese rind and pulp, which are characterized by different chemical properties.When dealing with heterogeneous samples, classical NIR spectroscopy may lead to inaccurate results since it is based on the acquisition of "average" spectra over a given sample area, thus losing the information related to the compositional variability within the sample.The importance of spatial information for the analysis of P-R grated cheese was recently demonstrated by Alinovi et al. (2019), which found a relationship between rind percentage and parameters related to particle size and distribution calculated from digital RGB images of the cheese samples.
The advantages of image-based methods and NIR spectroscopy can be coupled in NIR Hyperspectral Imaging (NIR-HSI), an analytical technique based on the acquisition of particular types of images, called hyperspectral images, where a whole NIR spectrum is registered for each image pixel (Gowen et al., 2007;Amigo et al., 2013;Calvini et al.;2018).More in detail, a hyperspectral image, also called hypercube, is a three-dimensional data array with two spatial dimensions (x pixel rows and y pixel columns) and one spectral dimension, corresponding to the λ acquired wavelengths.Therefore, the hypercube can be seen as a stack of spectrally resolved images, each one acquired at a given wavelength, or as a series of spatially resolved spectra, where each spectrum characterizes one pixel of the image (Wu and Sun, 2013a;Amigo et al., 2015;Calvini et al.;2016).
Considering that each NIR spectrum acts like a fingerprint of the chemical properties of a specific pixel, thanks to hyperspectral imaging it is possible to obtain both spatial and spectral information from a sample by observing variations in the chemical composition over the sample surface.These aspects resulted to be particularly useful for the analysis of heterogeneous matrices like food products (Wu and Sun, 2013b;Dale et al., 2013;Liu et al.,2017), generally overcoming the performances obtained with single-point spectroscopy (Burger and Geladi, 2006;Shonbichler et al., 2013).
However, the high amount of information contained in hyperspectral images can also become a drawback, since each hypercube can contain up to tens of thousands of pixel spectra, resulting in data handling and data storage issues.Indeed, in order to extract the relevant information from this kind of data, the application of proper multivariate statistical methods is mandatory (Burger and Gowen, 2011).This is known as Multivariate Image Analysis (MIA), which is based on the application to images of common chemometric methods, like e.g.Principal Component Analysis (PCA).This approach essentially consists in considering each pixel of the image as a separate object and the main goal is to find similarities or differences among clusters of pixels based on their spectral signatures (Prats-Montalban et al., 2011;Amigo et al., 2015).
When dealing with a large number of hyperspectral images that should be analysed altogether, classical pixel-level MIA can become unfeasible due to the intensive computational loads, since it would imply the simultaneous analysis of numerous images, each one with tens of thousands of pixel spectra.In these situations, a possible solution is to move from a pixel-level approach to an image-level approach, which consists in performing the analysis considering the image of each sample as a single object and extracting a feature vector characterizing the whole image, and thus the corresponding sample.In this manner, it is possible to analyse data matrices containing these feature vectors, in order to gain a general overview of the image dataset, to identify images sharing similar features or to quantify whole sample properties (Ulrici et al., 2012;Kucheryavskiy, 2013;Giraudo et al., 2018;Orlandi et al., 2018a;Orlandi et al., 2018b;Oliveri et al., 2019).
To this aim, a data dimensionality reduction method has been proposed, which consists in converting each hyperspectral image of the dataset into a one-dimensional signal, named hyperspectrogram, obtained by merging in sequence the frequency distribution curves of quantities derived from a PCA model calculated on the images (Ferrari et al., 2013;Ferrari et al., 2015;Xu et al., 2016;Calvini et al., 2016;Corti et al., 2017).In this manner, each hyperspectrogram summarizes the relevant information contained in the corresponding hyperspectral image and a large dataset of hyperspectral images is converted into a matrix of signals, which in turn can be analysed by means of common chemometric methods.
In this context, the main goal of the present study consisted in evaluating the possibility of exploiting the advantages of NIR-HSI coupled to data dimensionality reduction, in order to monitor the rind percentage of grated P-R cheese products.In particular, hyperspectral images of grated P-R cheese samples were analysed by means of the hyperspectrograms approach, and the resulting hyperspectrograms were then used to predict the rind percentage by means of Partial Least Squares (PLS) regression.
In order to minimize possible effects of unwanted variations, the mixtures were prepared starting from the same matrices of cheese pulp and rind, obtained by grating pulp and rind pieces derived from different cheese wheels matured for a period of 12 months.The mixtures were replicated twice (deliveries A and B) as reported in Table 1, each time following a different random order.Firstly, the matrices of grated pulp and rind were prepared, and a part of them was then used to obtain the first set of 15 mixtures.The remaining part of the grated pulp and rind matrices was stored in the dark at 4 °C, and after one week it was used to prepare the second set of 15 mixtures.For both the replicate sets, the samples were stored in the dark at 4°C and the day after their preparation they were delivered to the laboratory, where they were immediately analysed.
Furthermore, 15 additional samples with unknown rind percentage were provided by Parmigiano Reggiano cheese Consortium (X1-X15).These samples were prepared, delivered and stored considering the same procedure followed for the samples with known RP values.

Image acquisition
For each sample, three randomly sampled aliquots containing about 13 g of grated cheese were collected and placed inside a plastic Petri dish of 5.5 cm diameter.Each hyperspectral image included the three aliquots of two different samples.More in detail, the samples were positioned according to a 3  2 chessboard scheme, as reported in Figure 1.
The hyperspectral images were acquired using a line scanning system (NIR Spectral Scanner, DV Optic) equipped with a Specim Imspector N17E imaging spectrometer coupled to a Xenics Xeva-1.7-320camera (320  256 pixels) embedding Specim Oles 31 f/2.0 optical lens.The hyperspectral system covers the 900-1700 nm spectral range with a spectral resolution equal to 5 nm.Due to the low signal-to-noise ratio at the edges of the spectral range, only 143 spectral channels between 960 and 1670 nm were considered for further analysis.
The hyperspectral images were acquired using a black silicon carbide sandpaper sheet as background and including in the image scene also a white ceramic tile with a 99 % reflectance standard reference and two ceramic tiles with intermediate reflectance values corresponding to 89 % and 46 %, respectively.The raw data were converted into reflectance values using the instrumental calibration based on the measure of the white high reflectance standard reference and of the dark current (Ulrici et al., 2013).

Image elaboration
The first step of image elaboration consisted in the application of an additional internal calibration procedure in order to reduce possible variations over time.This correction procedure is based on the comparison of the average reflectance values of the white standard reference, of the two ceramic tiles and of the black silicon carbide sandpaper between an image chosen as reference and all the remainder images of the dataset.Further details about the image correction algorithm can be found in Ulrici et al. (2013).
Subsequently, the corrected images were cropped in order to obtain a single hyperspectral image for each aliquot of grated cheese sample.At the end of the cropping procedure, a total of 135 hyperspectral images were obtained (= 45 grated cheese samples  3 aliquots), as reported in the last column of Table 1.
After cropping, the pixels related to the black sandpaper background were removed using a thresholding procedure carried out considering a wavelength equal to 1100 nm.Indeed, at 1100 nm the pixels with reflectance value lower than 0.5 were ascribable to the background or to the plastic Petri dish, thus they were not considered in the subsequent steps.
Finally, an additional morphological erosion procedure was performed using a disk-shaped structuring element with radius equal to 8 pixels (Van Den Boomgaard and Van Balen, 1992).
Morphological erosion allowed to remove the pixels placed at the edges of the region of interest obtained after background removal, since these pixels were mainly influenced by scattering phenomena and specular reflections of the plastic Petri dish.
The image elaboration steps, which are summarized in Step 1 and Step 2 of Figure 1, were performed with routines written ad hoc in MATLAB language (ver.9.3, The Mathworks Inc., USA).

Exploratory analysis
As a preliminary assessment, Principal Component Analysis (PCA) was performed at the pixel level on some representative images (Prats-Montalbán et al., 2011).More in detail, three images corresponding to RP values equal to 0%, 20% and 40% were merged together and analysed by means of PCA after applying standard normal variate (SNV) and mean center as spectral preprocessing methods.This preliminary analysis was carried out in order to obtain a qualitative evaluation of the differences between samples containing an increasing amount of rind.

Data organization
Before calculating the calibration model to predict the rind percentage, the hyperspectral images of the grated cheese samples were split into training images, for model computation, and test images for external validation.The training images included the hyperspectral images of grated cheese samples with RP values equal to 0%, 10%, 14%, 18%, 22%, 26%, 30% and 40%, for a total of 48 images (= 8 RP values  2 deliveries  3 aliquots).The remainder images were separated into two different test sets: -TS known : including the images acquired on cheese samples with RP values equal to 5%, 12%, 16%, 20%, 24%, 28% and 35%, for a total of 42 images (= 7 RP values  2 deliveries  3 aliquots); -TS unknown : including the images acquired on the cheese samples of unknown composition, for a total of 45 images (= 15 unknown samples  3 aliquots).

Conversion into Common Space Hyperspectrograms
The hyperspectral images were then converted into one-dimensional signals, named Common Space Hyperspectrograms (CSH) (Calvini et al., 2016).The basic idea behind the hyperspectrograms approach consists in converting each hyperspectral image of the dataset into a one-dimensional signal, which acts like a feature vector retaining the useful spectral/spatial information of the corresponding image (Ferrari et al., 2013;Ferrari et al., 2015).More in detail, in the case of CSH, The first step in the computation of the CSH consisted in unfolding all the three-dimensional hyperspectral images into two-dimensional matrices with as many rows as the pixels retained after background removal and erosion, and as many columns as the number of spectral channels.Then, the unfolded hypercubes were row-preprocessed using SNV and scaled according to the global mean spectrum, obtained by averaging all the retained pixel spectra of the training images.After unfolding and spectra preprocessing, for each training image the corresponding variance-covariance matrix was calculated.Then, all the resulting variance-covariance matrices were summed in order to obtain the kernel variance-covariance matrix of the whole training set (Geladi and Grahn, 1996).
The kernel-variance covariance matrix was then decomposed by singular value decomposition (SVD) to obtain the loading vectors of the common PC space.In this case, the common PC space was calculated considering 3 principal components, based on the results of the previous exploratory data analysis described in Section 2.4.1.
Once calculated the PC space common to the training images, all the hyperspectral images belonging to both training and test sets were projected onto the PC space to obtain the corresponding scores, Q residuals and Hotelling's T 2 vectors.Finally, for each image, the corresponding CSH signal was obtained by merging in sequence the frequency distribution curves of the three score vectors, of the Q residuals vector and of the Hotelling's T 2 vector.
The frequency distribution curves were calculated considering a number of bins equal to 150 and normalized according to the number of pixels retained after image elaboration, as described in Section 2.3, to give the corresponding hyperspectrogram.Therefore, in this case each hyperspectrogram was a 750 points long vector, resulting from 150 bins 5 quantities derived from PCA (3 PCs + Q residuals + Hotelling's T 2 ).Further details about the algorithm used to calculate the CSH can be found in Calvini et al. (2016)

Calibration model
The training set matrix containing the CSH signals calculated from the training images was used to calculate the calibration model to predict the RP value, using Partial Least Squares (PLS) algorithm (Geladi and Kowalski, 1986).The signals were preprocessed using mean center and the optimal number of latent variables (LVs) was chosen by minimizing the cross-validation error.In particular, a custom cross-validation scheme was followed, considering 2 deletion groups, each one containing the signals derived from samples belonging to the same delivery day.
The performances of the calibration models were evaluated both in terms of Root-Mean-Square Error (RMSE) and of coefficient of determination (R 2 ).These parameters were calculated in calibration (RMSEC and R 2 Cal), cross-validation (RMSECV and R 2 CV) and prediction of the external test set (RMSEP and R 2 Pred).
The conversion of the hyperspectral images in CSH signals was done using a specific Graphical User Interface (GUI), that was previously developed by some of the authors of the present work.
The GUI, which works under the MATLAB environment (ver.9.

PCA at the pixel-level
For a first evaluation of the differences between grated cheese samples containing different amounts of rind, three hyperspectral images corresponding to samples with RP values equal to 0%, 20% and 40% were merged together to obtain a unique hyperspectral image, which was analysed at the pixellevel by PCA. Figure 3.a reports the resulting PC1-PC2 score plot, where each object represents a single pixel and is coloured according to pixel density, i.e. a yellowish colour represents a region of the PC1-PC2 score space with a high density of pixels, while blue corresponds to low pixel density.
From this score plot it is possible to observe the presence of three clusters of pixels, corresponding to the imaged samples with different RP values.The separation between samples with different rind levels is particularly evident along PC2.Indeed, the sample containing only cheese pulp is characterized by higher PC2 score values, while samples with increasing percentages of rind have decreasing PC2 score values, as shown in the PC2 score image reported in Figure 3.b.
In order to investigate the spectral features involved in the definition of the PC space, the corresponding PC1-PC2 loading vectors are reported in Figure 3.c.The highest absolute values of the PC2 loading vector can be found in the 1195-1225 nm wavelength range, corresponding to the C-H bond second overtone ascribable to lipids (Burns and Ciurzak, 2008;Karoui et al., 2006), in the 1330-1340 nm spectral range corresponding to asymmetric stretching vibration of water (Ozaky, 2002), and in the region centred at 1400 nm ascribable to the O-H bond first overtone of free water (Burns and Ciurzak, 2008).
Therefore, the amount of rind of the cheese samples can be somehow described by the distribution of the corresponding pixel spectra along the principal components.A convenient way to summarize this pixel distribution consists in using the frequency distribution curves of the score vectors of each sample, as reported in Figure 3.a that shows the frequency distribution curves of both PC1 and PC2 score vectors for each image.From this figure it is possible to observe the presence of a clear shift of the frequency distribution curves of PC2, which is related to the rind percentage of the corresponding samples.Although less marked, a variation with RP can be observed also for the frequency distribution curves of PC1 scores, which tend to become sharper and with a maximum located at lower PC1 values with increasing values of rind percentage.
Since the hyperspectrograms approach is based on the use of frequency distribution curves of score vectors calculated from a PCA model in order to summarize the relevant information contained in the images, this preliminary analysis suggests the effectiveness of this approach for the determination of the rind amount in hyperspectral images of grated cheese samples.

PLS calibration model with the CSH approach
The training set of CSH signals reported in Figure 2 was then used to calculate a PLS regression model for the quantification of rind percentage in the samples of grated Parmigiano Reggiano cheese, leading to the results reported in Table 2.
The optimal model dimensionality was found to be equal to 8 LVs, leading to a RMSECV value equal to 1.70 RP units, corresponding to a R 2 CV value of 0.979.
The calibration model was then used to predict the RP of the samples belonging to the TS known test set.In this case, the prediction results were calculated considering both the whole range of rind levels (0% -40%) and only the interval of rind percentages ranging from 10% to 30%, which better reflects RP values that generally may occur in real situations.
The prediction results confirm the effectiveness of the calibration model in quantifying the rind percentage, leading to RMSEP values equal to 1.91 and 1.85 RP units in the 0-40% and 10-30% ranges, respectively.3c.

Conclusions
The present study demonstrated the possibility of using NIR-HSI as a tool for the quantification of the amount of rind in grated Parmigiano Reggiano cheese samples.The combined use of a data dimensionality reduction approach, namely the Common Space Hyperspectrograms approach, with PLS regression allowed to obtain satisfactory prediction performances.

Conflict of interest
The authors declare that they have no conflict of interest.

*Conflict of Interest Form
the signals are obtained by merging in sequence the frequency distribution curves of quantities derived from a common PCA model, i.e. from a model calculated considering all the images of the training set.
. The conversion of the hyperspectral images into CSH signals is schematically depicted in Step 3 of Figure 1.In this manner, at the end of the conversion procedure three different matrices of signals were obtained: the training set (TR), the test set derived from the TS known images and the test set of unknown samples derived from the TS unknown images.

Figure 2
Figure 2 shows a plot of the CSH signals belonging to the training set, coloured according to the rind percentage of the corresponding sample.

Figure 4
Figure4shows the plot of the rind percentage values predicted for the TS known test set versus the corresponding experimental values, where the samples are coloured according to the delivery day.Generally, all the objects are close to the bisector, indicating the good prediction performances of the model.In addition, from Figure4it is also possible to observe that there are not evident systematic variations in the prediction results caused by the different delivery days, which further confirms the robustness of the calibration model toward replicated measurements.Considering that compliant Parmigiano Reggiano grated cheese samples should have an RP value lower or equal than 18%, all the test set samples with a rind percentage falling outside this limit are correctly identified by the model.The calibration model was also used to predict the samples with unknown composition belonging to the TS unknown test set.The predicted rind percentages of the unknown samples were communicated to the Parmigiano Reggiano Cheese Consortium, which then revealed the corresponding experimental values.In this manner, it was possible to perform a further external validation of the calibration model.It has to be considered that, for each unknown sample, three different aliquots were imaged and, therefore, three RP values were obtained in prediction by the model.In order to have a single estimate of the RP value for each unknown sample, the three RP predicted values corresponding to the three aliquots were averaged for each sample.The results are reported in Table3, together with

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Key steps involved in image elaboration and analysis

Figure 4 .
Figure 4. Results of the PLS model: TSknown test set predicted rind percentage (Y Predicted) vs experimental rind percentage (Y measured).

Figure 5 .Highlights
Figure 5. VIP scores of the PLS model (a) and PC2 score images of samples with increasing rind percentage values(b).

Figure 5
Figure 5 Table3, together with the corresponding experimental values, further confirming the good prediction ability of the calibration model.More in detail, an RMSEP value equal to 2.50 RP units was obtained, corresponding to an R 2 value equal to 0.955.The highest difference between predicted and experimental rind percentage, equal to 5 RP units, was observed for sample X15, while samples X2 and X4 were exactly predicted.In order to have a comprehensive evaluation of the hyperspectrogram regions most relevant to the calibration model, Figure5.areports the Variable Importance in Projection (VIP) scores: variables with VIP score values higher than one are considered significant for the model.This figure shows that all the frequency distribution curves of the PCA quantities included in the CSH signals have regions with significant variables.In particular, among the three score vectors included in the signals, the regions related to the frequency distribution curve of PC2 reach the highest VIP score values, together with the signal regions related to Hotelling T 2 vales.As an example, Figure 5.b reports the PC2 score images of some representative test set samples with increasing percentage of rind.Similarly to what was previously reported in Section 3.1, images of samples with a lower RP value are characterized by higher PC2 score values, and the increase of the RP value in the grated cheese samples causes a shift of the corresponding pixels toward lower PC2 score values.Actually, this is due to the fact that the PC2 loading vector of the common PCA model used to calculate the CSH signals is very similar to the PC2 loading vector calculated on the three sample hyperspectral images, reported in Figure

Table 1 .
The calibration model was validated using two different test sets: the first test set consisted of cheese samples with known RP values, while for the second test set the experimental RP values of the analysed samples were revealed by the operators of Parmigiano Reggiano Cheese Consortium only after providing them with the RP values predicted by the model.The RMSEP values obtained for both test sets, Summary information about the grated cheese samples considered in the present study processing methods, among others.Therefore, in order to obtain a more robust calibration model the influence of these factors should be properly evaluated and included in the model.As aCaptions to Tables and Figures

Table 2 .
Results of PLS regression for the determination of rind percentage.

Table 3 .
Prediction results of the unknown test samples and corresponding experimental RP values.

Table 1 .
Summary information about the grated cheese samples considered in the present study

Table 2 .
Results of PLS regression for the determination of rind percentage.

Table 2 Table 3 .
Prediction results of the unknown test samples and corresponding experimental RP values.