## 2 Materials and Methods
### Data Analysis
#### 2.4.2 CCP Analysis  
A 3-dimensional DCNN (3D-DCNN) model was implemented to classify the thermal images recorded during the whole CCP (i.e., baseline and recovery conditions) in the three classes (HC, SSc, and PRP). Due to the high computational load and the time required to analyze all the images of the recorded IR video (one every 30 s for 20 min, 40 images in total), a representative sample was considered for the successive analyses. This sample was composed of IR-images extracted every three minutes. The baseline image was also added to the set of images to better describe the rewarming phenomenon and in accordance with previous studies [13; 19]. Figure 4: Group average values and the related standard deviation of each feature z-score normalized. The measurement units reported on the y-axis are the standardize units since the features shown are z-score normalized. The 3D-DCNN differs from the 2D in the convolutional filter type. Indeed, it applies a three-dimensional filter to the dataset and the filter moves in three-directions (x, y, z) to calculate the low-level feature representations. Their output shape is a three-dimensional volume space. Architectures with volumetric (i.e., spatially 3D) convolutions have been successfully used in video analysis or 3D medical images [45; 46]. In this case, time acts as the third dimension [47]. The 3D-DCNN model was implemented following the same processing step used to develop the DCNN model (Section 2.4.1). The batch of input images underwent the same preprocessing procedure as described in Section 2.4.1, whereas the model architecture design slightly changed. The number of 3D-filters in the three convolutional layers remained the same, whereas the filter sizes were, respectively, 5 \(\times\) 5 \(\times\) 4, 3 \(\times\) 3 \(\times\) 2, 3 \(\times\) 4 \(\times\) 3, and the three MaxPooling, respectively, 3 \(\times\) 3 \(\times\) 1, 3 \(\times\) 3 \(\times\) 1, 3 \(\times\) 4 \(\times\) 2. In this way the dimensionality of each feature was reduced to 1 before the fully connected layer. Finally, the model optimization and evaluation were conducted as previously described. The features map resulting as output from the first 2 convolutional layers together with the batch of input images, respectively for HC, SSc, and PRP participants are shown in Figure 5.  
Figure 5: 3-Dimensional Deep Convolutional Neural Network (3D-DCNN) feature maps. (**a**) Batch of input images representative of the cold challenges procedure (CCP) data of a healthy participant randomly chosen, together with the features map resulting from the first two 3D-convolutional layers. (**b**) Batch of input images of a systemic sclerosis (SSC) participant randomly chosen, and the related features map resulting from the first two 3D-convolutional layers. (**c**) Batch of input images of a primary Raynaud’s phenomenon (PRP) participant randomly chosen, and the related features map resulting from the first two 3D-convolutional layers. The feature-based analysis on the CCP data was performed by extracting features of interest from the fingers' temperature trend. Such a trend was obtained by averaging the temperatures of the selected regions of interest (see Figure 2) throughout the experiment. The temperature time courses of each finger of a representative HC, SSc, and PRP participant are shown in Figure 6.  
The extracted temperature data were then normalized by subtracting the average value of the entire dataset. The features identified to characterize the dynamic response of the finger's temperature to the cold challenge were the following:  
1. Delay time: the time required for the finger's temperature after the cold challenge to reach 50% of its final value (i.e., the recovery temperature after 20 min from the cold stress). 2. Rise time: the time required for the finger's temperature to rise from 10% to 90% of its final value. 3.