Automatic marbling prediction of sliced dry-cured ham using image segmentation, texture analysis and regression

Expert Systems With Applications 206 (2022) 117765 Available online 11 June 2022 0957-4174/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Expert Systems With Applications journal homepage: www.elsevier.com/locate/eswa Automatic marbling prediction of sliced dry-cured ham using image segmentation, texture analysis and regression Eva Cernadasa,∗, Manuel Fernández-Delgadoa, Elena Fulladosab, Israel Muñozb aCentro Singular de Investigación en Tecnoloxías Intelixentes da USC (CiTIUS), Universidade de Santiago de Compostela, Rúa Xenaro de la Fuente Domínguez, Santiago de Compostela, 15782, Spain bIRTA-TA, Food processing and Engineering, Finca Camps i Armet, Monells, Girona, 17121, Spain ARTICLE INFO Keywords: Dry-cured ham Intramuscular fat Marbling Support vector regression Texture analysis Image segmentation ABSTRACT Dry-cured ham is a traditional Mediterranean meat product consumed throughout the world. This product is very variable in terms of composition and quality. Consumer’s acceptability of this product is influenced by different factors, in particular, visual intramuscular fat and its distribution across the slice, also known as marbling. On-line marbling assessment is of great interest for the industry for classification purposes. However, until now this assessment has been traditionally carried out by panels of experts and this methodology cannot be implement in industry. We propose a complete automatic system to predict marbling degree of dry-cured ham slices, which combines: (1) the color texture features of regions of interest (ROIs) extracted automatically for each muscle; and (2) machine learning models to predict the marbling. For the ROIs extraction algorithm more than the 90% of pixels of the ROI fall into the true muscle. The proposed system achieves a correlation of 0.92 using the support vector regression and a set of color texture features including statistics of each channel of RGB color image and Haralick’s coefficients of its gray-level version. The mean absolute error was 0.46, which is lower than the standard desviation (0.5) of the marbling scores evaluated by experts. This high accuracy in the marbling prediction for sliced dry-cured ham would allow to deploy its application in the dry-cured ham industry. 1. Introduction Dry-cured ham is a traditional meat product of many Mediterranean countries that is widely consumed throughout the world, being its flavor and texture characteristics highly appreciated by consumers. There are many factors affecting the final characteristics of dry-cured ham, such as processing conditions and raw material characteristics, i.e. fat content (Coll-Brasas et al.,2021). In sliced dry-cured ham, visual intramuscular fat, subcutaneous fat thickness and color are the parameters most used by the consumers for the product evaluation, therefore affecting consumer’s acceptability and purchase decision (Lorido et al., 2019). Although several non-invasive technologies can be used to categorize entire hams according to its fat content (de Prados et al., 2015), these technologies cannot be used to predict intramuscular fat (IMF) in sliced products because of its variability between the muscles of a ham. Eating quality in meat has been associated to the fat distribution rather than to the total IMF. Distribution of IMF is usually known as marbling, because of its appearance similar to marble (Cernadas ∗Corresponding author. E-mail addresses: [email protected] (E. Cernadas), [email protected] (M. Fernández-Delgado), [email protected] (E. Fulladosa), [email protected] (I. Muñoz). et al.,2002). One of the most important challenges for producers is the heterogeneity of the marbling in slices, that can vary significantly among ham pieces and even within the same piece. Classification of slices of dry-cured ham according to the marbling degree is of special interest for the food industry. Producers would be able to segment the market, offering products tailored to consumer’s needs and increasing the value of their production. Marbling ranking in different meats and meat products has been performed by panels of trained experts or relying on standards consisting of pictures depicting a scale of marbling (from 0.0 to 10.0), as it is the case for the National Pork Producers standards (Moines,1999). In the case of dry-cured ham, a marbling ranking has been developed but it is not still published. However, marbling evaluations by experts are costly and are not feasible for the ham industry. Computer image analysis might be a solution because it is a fast and non-destructive technology, and it is a replicable and repetitive method that has been successfully applied to the assessment of multiple food characteristics: fish (Dutta et al.,2016), cheese (Dias et al.,2021) or bread (Srivastava https://doi.org/10.1016/j.eswa.2022.117765 Received 26 October 2021; Received in revised form 2 May 2022; Accepted 3 June 2022 Expert Systems With Applications 206 (2022) 117765 2 E. Cernadas et al. et al.,2015). The scientific literature includes several studies that apply computer image analysis to determine IMF and marbling in different meat products. Combining magnetic resonance imaging (MRI) and computer vision techniques, the works (Ávila et al.,2019;Cernadas et al.,2005) predicted marbling in the biceps muscle of dry-cured hams and loins, but MRI is a technology that is expensive to install in meat industries. For the segmentation of IMF in meat and meat products, several techniques have been applied: K-means clustering in beef Longissimus dorsi muscle (Jackman et al.,2009), automatic thresholding (Liu et al.,2018) and the Otsu method (Uttaro et al.,2021) in pork loin, multi-scale line detection (Cernadas et al.,2002), gradient-based techniques (Santos-Garcés et al.,2014) and convolutional neural networks (Muñoz et al.,2019) for dry cured ham. In general, these segmentation algorithms perform quite well, high correlation values or low classification errors are obtained, depending on the aim of the study. For classification of marbling, the following techniques have been applied: chemical pre-treatments and line detection algorithms (Faucitano et al.,2004), line detection algorithm (Huang et al.,2013;Liu et al.,2012) in pork meat, neural networks (Muñoz et al.,2015) and decision trees in hyperspectral images of dry cured ham (Velásquez et al.,2017). In general, the results of marbling classification are quite satisfactory with low prediction errors and at least 90% of the samples correctly classified. However, the evaluation of IMF and marbling in dry-cured ham slices is still a challenge. A wide range of variation in the color of lean and fat tissues can be observed in slices across hams, which poses a challenge for image segmentation and the evaluation of marbling. These differences are explained by the different levels of drying of the hams, and the presence of precipitates such as phosphates of tyrosine crystals, with a white color similar to that of the fat. This paper proposes a prototype to automatically predict marbling of the principal muscles from a ham slice using image segmentation, texture analysis and regression models. Specifically, we define an algorithm that automatically extracts squared regions inside the main muscles. Then, color texture features are computed for each region, which are the inputs to a regression model that predicts the marbling score for each muscle. The paper is organized as follows. Section 2 describes the materials used to obtain the ham slice, to develop the sensorial analysis on the ham muscles and to annotate the contour of ham muscles. Section 3describes the algorithm used to extract automatically the square ROIs representing each ham muscle, and briefly explains the color texture features extraction techniques and the regression models. Section 4describes the experimental setup and statistical evaluation measures used. Section 5presents and discusses the results. Finally, Section 6summarizes the main conclusions and proposals of future work. 2. Materials This section describes the material used to obtain the data in this research: the system used to capture images of ham slice (Section 2.1), the traditional procedures to determine the marbling of a ham muscle (Section 2.2), and the process to draw the outline of each muscle on the ham slice (Section 2.3). 2.1. Image acquisition High quality images were acquired with a calibrated digital camera Canon EOS 50D (15.1 megapixels) and an objective Canon EF-S 18– 200 mm f/3,5–5,6 IS. The camera was mounted in a black closet (1.06 ×1.06 ×2.50 m3) with 8 equidistant halogen lights Solux Q50MR16 CG/47/36◦12 V/50 W/4700 K (Eiko Ltd., Shawnee, Kansas, U.S.A.) to ensure a correct lighting. White balance was carried out with a white card (Lastolite). The camera was connected to a computer to store the images. Slices were placed 30 cm below the camera on a uniform black Fig. 1. Scheme of the image acquisition system. Table 1 Number of images, minimum, maximum, average and standard desviation of the marbling values for the different ham muscles used in this experimentation. Muscle #images Min. Max. Avg. Dev. Biceps femoris 337 1 7 3.1 1.0 Semimembranosus 322 0.5 6 2.0 1.0 Semitendinosus 55 4 9 6.2 1.2 surface. Photos of both sides of the dry-cured ham slices were taken. All the images were taken during the same session. The white balance of the images was carried out with Capture One PRO 5.0 software (Phase One A/S Inc., Frederiksberg, Denmark) and RGB images of 667 ×1000 pixels with 16 bits color were obtained, with one pixel corresponding to 0.3968 mm2. For the evaluation of marbling, the computer screen was calibrated so that the colors of the images were as close as possible to the colors of the samples (NEC Multisync LCD 2690 WUXI2). Fig. 1 shows a scheme of the image acquisition system used. 2.2. Marbling evaluation Sensory analysis (marbling evaluation) of the samples was carried out by six trained panelists (ISO 8586-2: 2012) and consisted of a visual assessment of the marbling of the most representative muscles: Biceps femoris (BF), Semimembranosus (SM) and Semitendinosus (ST) of a dry-cured ham (Bermúdez et al.,2014). Marbling was scored by consensus (in our case three panelists) by means of scoring scale from 0.5 (minimum marbling) to 10 (maximum marbling) at intervals of 0.5. When scoring marbling, the panelists considered the total amount and the distribution of the fat streaks. Marbling evaluation was done in triplicate by the panelists. The standard deviation of the panelists among trials was determined at 0.5 points. A collection of commercial dry-cured hams were obtained from different ham producers with crosses from different pig breeds (Large White, Landrace, Duroc and Iberian) and having a wide range of marbling. A 2 cm thick slice containing muscles BF, SM and ST was obtained at 10 cm from the aitch bone in the distal direction (at the widest part of the ham) and packed into plastic bags of polyamide/polyethylene (oxygen permeability of 50 cm3/m2/24 h at 23 ◦C and water permeability of 2.6 g/m2/24 h at 23 ◦C and 85% RH, Sacoliva©S.L., Spain). The image dataset is composed of 714 images obtained in the following way: photos were obtained from 180 commercial dry-cured hams, 2 slices/ham (at different points in the ham, obtaining slices quite different one from another) and 2 muscles for each slice, giving a total of 180 hams ×2 slices ×2 muscles =720 photos. Six of these photos were not included in the evaluation due to defects on the surface such as cuts and phosphate crystals. For each image, it was only provided the measure of marbling for one muscle with values between 0.5 and 9with the distribution shown in Table 1. Fig. 2 shows examples of different marbling scores for biceps muscle. Expert Systems With Applications 206 (2022) 117765 3 E. Cernadas et al. Fig. 2. Examples of marbling scores for biceps muscle: (a) slice 8372 with score 1; (b) slice 8485 with score 1; (c) slice 8274 with score 3; and (d) slice 8424 with score 6.5. Fig. 3. Examples of the contours of biceps (left panel) and semimembranosus (right panel) muscles overlapped to images of ham slice. Fig. 4. Stages of the method to predict the marbling from dry-cured ham slices. 2.3. Muscle annotation To develop the first experiment, the experts delineated the contour of the muscle for which the marbling was estimated by sensorial analysis. To draw the contours for all images, they used a home-made software, as it can be seen in Fig. 3. 3. Methods The system proposed to predict the marbling from dry-cured ham slices, shown in Fig. 4, encloses the following stages: (1) the image acquisition system already described in Section 2.1; (2) the automatic extraction of the ROIs in the ham slice; (3) the computation of features from the ROI extracted; and (4) the regression model to predict the marbling score of each muscle in the ham slice. The Sections 3.1,3.2 and 3.3 describe the stages 2, 3 and 4, respectively. 3.1. Automatic extraction of ROIs The ham slice images are processed to automatically extract square ROIs from the biceps femoris, semimembranosus and semitendinosus muscles. These extracted ROIs will be used in the third experiment. In this process, we take into account the anatomical information about the distribution of the muscles and subcutaneal/intermuscular fat within the ham. As it can be seen in Fig. 3, some slices present a hole in the slice (left panel), due to the slice is cut by the ham bone, and others not (right panel). As well, the biceps muscle can be in the right or left side of the image. The algorithm to extract the square ROIs encloses the following steps: (1) extract the ham slice from the image; (2) check if the slice has the bone hole; (3) if there is not a hole in the slice, find the biggest intermuscular fat region in the slice; (4) in both previous cases (step 2 or 3), a reference position is calculated to know the slice orientation, which allows to know if the muscles are upper/bottom or if the BF muscle is on left/right side of the slice; and (5) extract a square region for each muscle. In our case, we use a ROI with length 𝑠= 64 pixels. To extract the ham slice, the original RGB image 𝐼(𝑥, 𝑦), with 𝑥= 1,…, 𝑁, and 𝑦= 1,…, 𝑀, of dimensions 𝑁×𝑀, is transformed to the Lab color space, because it is more robust to illuminance variance than the RGB space (Cernadas et al.,2017). Let 𝐼𝑏(𝑥, 𝑦)be the 𝑏channel of the ham slice after smoothing with a mask (we use a mask of 5 pixels) in order to attenuate the small fat features and noise (column b in Fig. 5). Let ℎ𝑏be the histogram of image 𝐼𝑏(𝑥, 𝑦). The maximum value 𝐻𝑏of ℎ𝑏is chosen to calculate the area of ham slice. The 𝐼𝑏(𝑥, 𝑦)image is Expert Systems With Applications 206 (2022) 117765 4 E. Cernadas et al. Fig. 5. Examples of the extraction of square ROIs from the ham slice for all muscles: (a) original images; (b) channel 𝑏of Lab color space after smoothing; (c) binary mask with the ham slice and bone hole and (d) the extracted ROIs for each muscle, white for biceps, yellow for semimembranosus and green for semitendinosus. Fig. 6. Extraction of the square ROIs of each muscle. Upper left panel: scheme of a ham slice. Upper right panel: row 540 of the transformed version 𝐼𝑎of the ham slice no. 8607. Lower panel: the ham slice with the extracted ROIs overlapped: biceps (white), semimembranosus (yellow) and semitendinosus (green). thresholded to calculate the binary image 𝐵(𝑥, 𝑦)using the following expression: 𝐵(𝑥, 𝑦) = {0|𝐼𝑏(𝑥, 𝑦) − 𝐻𝑏|≤5 255 otherwise (1) After thresholding, we apply morphological operators to the binary image 𝐵in order to fill small holes. First, the biggest region of 𝐵is extracted and it is associated with the contour of ham slice. Next, the algorithm searches for a large black region inside this contour. If this region is found, it is associated to the bone hole. The contours of ham slice and bone hole are used to create a binary mask image 𝐵𝑚(𝑥, 𝑦) = 0 (black color in column c of Fig. 5) when (𝑥, 𝑦)is inside the contour and outside the bone, and 𝐵𝑚(𝑥, 𝑦) = 1 (white color) when (𝑥, 𝑦)is outside the contour or inside the bone. This process is shown in the columns a, b and c of Fig. 5. The process to extract the ROI for each muscle, denoted as 𝐼𝐵𝐹 ,𝐼𝑆𝑀 and 𝐼𝑆𝑇 for biceps femoris,semimembranosus and semitendinosus, respectively, is drawn in Fig. 6 and summarized by Expert Systems With Applications 206 (2022) 117765 5 E. Cernadas et al. Algorithm 1: Automatic extraction of square ROI images for each muscle from a ham slice. 1Algorithm: [𝐼𝐵𝐹 , 𝐼𝑆𝑀 , 𝐼𝑆𝑇 ]=ExtractSquaredROI(𝐼, 𝑠) Data: 𝐼: original RGB image of ham slice; 𝑠: size of ROI Result: 𝐼𝐵𝐹 , 𝐼𝑆𝑀 , 𝐼𝑆𝑇 : square ROI images for biceps, semimembranosus and semitendinosus muscles 2𝐼𝑏←𝑏channel of Lab image smoothed by box filter 3ℎ𝑏←histogram of 𝐼𝑏 4𝐻𝑏←maximum of ℎ𝑏 5𝐵←binary image using Eq. (1) and morphological processing 6𝐵𝑚←image mask with ham slice outline and hole if exist 7𝑅𝑠←(𝑥𝑠, 𝑦𝑠, 𝑤𝑅𝑠, ℎ𝑅𝑠)rectangle enclosing ham slice 8𝐼𝑎←𝑎channel of Lab image smoothed and masked by 𝐵𝑚 9𝜇𝑎←average value of 𝐼𝑎inside 𝐵𝑚;offset←10 10 if existsHole(𝐼𝑏)then 11 𝑅ℎ←(𝑥ℎ, 𝑦ℎ, 𝑤𝑅ℎ, ℎ𝑅ℎ)rectangle enclosing hole 12 𝑑𝑢𝑝𝑝𝑒𝑟 ←𝑦𝑠−𝑦ℎ;𝑑𝑙𝑜𝑤𝑒𝑟 ←𝑦𝑠+ℎ𝑅𝑠− (𝑦ℎ+ℎ𝑅ℎ) 13 if 𝑑𝑢𝑝𝑝𝑒𝑟 > 𝑑𝑙𝑜𝑤𝑒𝑟 then 14 𝑦𝑐←𝑦𝑠− (𝑦𝑠−𝑦ℎ)∕4 15 𝑦𝑒←𝑦ℎ−𝑠∕2 − offset 16 else 17 𝑦𝑐←𝑦𝑠+ℎ𝑅𝑠−𝑠− (𝑦𝑠+ℎ𝑅𝑠−𝑦ℎ−ℎ𝑅ℎ)∕4 18 𝑦𝑒←𝑦ℎ+ℎ𝑅ℎ+𝑠∕2 + offset 19 end 20 𝑥𝑐←middlePointX(row(𝑦𝑐+𝑠∕2)) 21 𝐼𝑆𝑇 ←extST(𝑥𝑐−𝑠∕2,𝑦𝑐,𝑠) 22 [𝐼𝐵𝐹 , 𝐼𝑆𝑀 ]←extBSM(𝑦𝑒) 23 else 24 𝑇2←second Otsu’s threshold of 𝐼𝑎 25 (𝑥𝑖𝑛, 𝑦𝑖𝑛)←centroid of the largest inner region inside 𝐼𝑎 after thresholding using 𝑇2 26 𝐼𝑆𝑇 ←extSTC(𝑥𝑖𝑛, 𝑦𝑖𝑛, 𝑠) 27 [𝐼𝐵𝐹 , 𝐼𝑆𝑀 ]←extBSM(𝑦𝑖𝑛) 28 end algorithm 1. Let 𝑅𝑠be the rectangle enclosing the ham slice, which is defined by the top left vertex (𝑥𝑠, 𝑦𝑠)and by its width (𝑤𝑅𝑠) and height (ℎ𝑅𝑠). Let 𝐼𝑎be the 𝑎channel of Lab image of the ham slice masked (multiplied) by 𝐵𝑚. Two cases can be considered: (a) When the ham slice has visible bone hole, the algorithm uses as reference its enclosing rectangle 𝑅ℎ, with initial coordinates (𝑥ℎ, 𝑦ℎ)and width and height 𝑤𝑅ℎand ℎ𝑅ℎ, respectively. (b) When the ham slice has not a visible bone hole, the algorithm searches for the largest intermuscular fat region and calculates its centroid. In order to discard the background from the image 𝐼𝑎, we selected a threshold applying the Otsu’s method (Otsu, 1979) with three thresholds {𝑇𝑖}3 𝑖=1, that correspond to different types of materials in the image (background, subcutaneal fat and muscle). Experimentally, we checked that threshold 𝑇2ensures an accurate segmentation of background from the remaining materials. In the binary image generated by thresholding 𝐼𝑎 with 𝑇2, the inner biggest region is selected as representing the intermuscular fat region, whose centroid is (𝑥𝑖𝑛, 𝑦𝑖𝑛). These reference points (bone hole or centroid of intermuscular fat) allow to select the height where the different muscles should be extracted. In order to extract the BF and SM muscles, the algorithm finds out whether the biceps is on the left or right side of the ham slice. The position in the horizontal axis for extracting the BF and SM muscles is determined analyzing a specific row in the image 𝐼𝑎. This row is smoothed in order to attenuate the random noise and it is denoted as 𝚛𝚘𝚠(𝑦𝑏) = 𝐼𝑎(𝑥, 𝑦𝑏),𝑥= 1,…, 𝑀 (see upper left panel of Fig. 6). In order to locate the BF muscle, we estimate the pixels representing the subcutaneal fat (close to the biceps muscle) along the 𝚛𝚘𝚠(𝑦𝑏)counting the number of values. Specifically, the procedure is as follows: 1. Let 𝑘1be the first value of 𝑥where 𝐼𝑎(𝑥, 𝑦𝑏)>0coming from left to right (see the upper right panel of Fig. 6). Let 𝑤𝑥1be the number of values of 𝑥where 𝐼𝑎(𝑥, 𝑦𝑏)< 𝜇𝑎, for 𝑥=𝑘1, 𝑘1+ 1,…, 𝑀 (i.e. going from left to right side), where 𝜇𝑎is the average value of 𝐼𝑎inside 𝐵𝑚. 2. Let 𝑘2be the first value of 𝑥where 𝐼𝑎(𝑥, 𝑦𝑏)>0coming from right to left, and 𝑤𝑥2the number of values that 𝐼𝑎(𝑥, 𝑦𝑏)< 𝜇𝑎, for 𝑥=𝑘2, 𝑘2− 1,…,2,1(i.e. going from right to left side). If 𝑤𝑥1> 𝑤𝑥2the BF muscle is on the left side of ham slice. Otherwise, it is on the right side. This process is performed by the extBSM(y) function in the algorithm 1, where 𝑦represents the row to be analyzed in the image. So the square ROIs for muscles BF and SM, of size 𝑠, are extracted at positions (𝑘1+𝑤𝑥1 + 𝑠∕2 + offset, 𝑦) and (𝑘2−𝑤𝑥2 − 𝑠∕2 − offset, 𝑦)(we use an offset of 10 pixels in order to avoid defects in the contour of ham slice). The regions extracted for each muscle are shown overlapped to the ham slice in the column dofFig. 5. When there is bone hole, the 𝑦coordinate for extracting muscles are determined calculating the distances 𝑑𝑢𝑝𝑝𝑒𝑟 =𝑦𝑠−𝑦ℎand 𝑑𝑙𝑜𝑤𝑒𝑟 =𝑦𝑠+ℎ𝑅𝑠− (𝑦ℎ+ℎ𝑅ℎ). If 𝑑𝑢𝑝𝑝𝑒𝑟 > 𝑑𝑙𝑜𝑤𝑒𝑟, the muscles are above the bone hole, otherwise the muscles are below the bone hole. The middlePointX(y) function returns the middle point of the ham slice for the image row 𝑦. The extST(x,y,s) function extracts a square ROI of size 𝑠for the ST muscle from the original RGB image in the point (𝑥, 𝑦)for ham slices with visible hole. The extSTC(x,y,s) function extracts a square ROI of the ST muscle for ham slices without visible hole. Let 𝑥𝑐1=middlePointX(row(𝑦𝑖𝑛 +𝑠∕2)) and 𝑥𝑐2= middlePointX(row(𝑦𝑖𝑛 −𝑠∕2)) be two middle points in the ham slice in the horizontal axis near the centroid of intermuscular fat. Let 𝑟1and 𝑟2be two ROIs of size 𝑠, extracted from the image 𝐼𝑎in the points (𝑥𝑐1−𝑠∕2, 𝑦𝑠−(𝑦𝑠−𝑦𝑖𝑛)∕6) and (𝑥𝑐2−𝑠∕2, 𝑦𝑠+𝑤𝑅𝑠−𝑠−(𝑦𝑠+𝑤𝑅𝑠−𝑦𝑖𝑛)∕3) respectively. To select which ROI corresponds to the ST muscle, the mean value of both ROIs, 𝜇𝑟1and 𝜇𝑟2, are calculated and the ST muscle is the ROI with the highest mean value, which corresponds with the ROI containing more fat. 3.2. Color texture features Texture analysis has been employed in previous works to predict different qualities or attributes of meat products from MRI images (Ávila et al.,2019;Cernadas et al.,2005), which are gray level images. Nevertheless, it is known that the color is also a very important characteristic in other computer vision applications (Cernadas et al.,2017; González-Rufino et al.,2013). Color texture analysis can be tackled from different paradigms: simple color features, gray level texture features and integrative color texture analysis. A recent work (Cernadas et al.,2017) concluded that parallel approaches, that concatenate the two former, are superior than analyzing directly the color texture with integrative approaches. There are many methods to extract only the color in a strict sense (Cernadas et al.,2017). In this work we use two color spaces, RGB and Lab, where the chromatic channels are aand b. We use first-order features of three types, each with 2 feature vectors: 1. Only mean value for each chromatic channel, denoted as CM, that stands for ‘‘color mean’’: (1) CMRGB, 3 features: mean color of the channels R, G, and B of the muscle; (2) CMab, 2 features: mean color of the channels aand bof the Lab image. In both cases, the mean values are only calculated inside the muscle regions. Expert Systems With Applications 206 (2022) 117765 6 E. Cernadas et al. 2. Mean and variance for each chromatic channel, CMV, or color mean and variance: (1) CMVRGB, 6 features: mean and variance of the channels R, G and B of the RGB image. (2) CMVab, 4 features: mean and variance of the channels aand bof the Lab image. 3. First order statistics, denoted as FOS, for each chromatic channel, that include mean, variance, skewness, kurtosis and entropy: (1) FOSRGB, 15 features: 5 features ×3 color channels of the RGB image. (2) FOSab, 10 features: 5 features ×2 color channels of the Lab image. The most popular gray level texture features belong to the families of statistical, such as second-order features and local binary patterns (LBP), and frequencial, including wavelet and Gabor features. The gray level version of the original image is obtained following two alternative approaches: (1) converting a RGB image to gray level image; and (2) using the L channel of a Lab image. Among the second-order statistics we selected the Haralick coefficients (Haralick et al.,1973), derived from the gray level coocurrence matrix (GLCM). These coefficients describe the probability of finding two pixels with the same value at different scales, or distances, and orientations, or angles. The parameters used normally are: (1) the orientations 0◦,45◦,90◦and 135◦; and (2) three scales, with pixel distances of 1, 2, and 3. For each scale, the GLCM matrix is averaged over all orientations, and the contrast, homogeneity, correlation and energy of the matrix is computed. Two feature vectors were considered: (1) HarRGB, that includes the four previous features for scales {1,2,3} calculated on the gray version of the RGB image, with 4 ×3=12 features; and (2) HarLab, which is similar to HarRGB but calculated on the L channel of Lab image. When the input is an irregular region, the features are computed only on pixels included in the region (González-Rufino et al.,2013). The Haralick’s coefficients were computed using the graycomatrix() function of the Matlab Image Processing Toolbox.1 The LBP operator is a state-of-art texture analysis approach proposed by Ojala et al. (2002), which describes each pixel comparing its value with the neighboring pixels. For each neighboring pixel, the result will be set to one if its value is higher than the value of central pixel, otherwise the result will be set to zero, developing a binary code for each pixel. We use the uniform LBP, which considers the binary patterns with only two transitions (from 0 to 1 and vice versa). In a circularly symmetric neighbor set of 𝑃pixels can occur 𝑃+ 1 uniform binary patterns. The number of ‘‘1’s’’ in the binary pattern is the label of the pattern, while the nonuniform patterns are labeled by 𝑃+ 1. This process can be applied to different scales, as in the case of coocurrence image. The histogram of the pattern labels accumulated over the intensity image is employed as texture feature vector. The most common values for these parameters are {(𝑃 , 𝑅) ∈ (8,1),(12,2),(16,3)}, where 𝑃is the number of neighbors and 𝑅is the distance between the central pixel and the neighbors. We construct the texture feature vectors LBPRGB and LBPLab to be applied on the gray version on RGB image and on the L channel of the Lab image, respectively. Both vectors have 42 = (8 + 2) + (12 + 2) + (16 + 2) features, because the uniform LBP are 𝑃+2 features. We use the LBPMatlab2implementation provided by the LBP creators. Discrete wavelet transform (DWT) representation is a theory for multi-dimensional signal decomposition (Laine & Fan,1993;Walker, 2008) which recursively apply filters to decompose the image into low-pass and high-pass frequency bands. A compact representation for texture analysis can be computed taking the mean and variance of the energy distribution for the transformed coefficients in each sub-band and decomposition level. We compute texture feature vectors calculating the mean and variance of the energy over 3 levels of decomposition 1https://www.mathworks.com. 2http://www.cse.oulu.fi/CMV/Downloads/LBPMatlab. and create the feature vectors 𝐖𝑇𝑖𝑗𝑘, where: (1) 𝑖= {𝐻𝑎𝑎𝑟, 𝐷𝑎𝑢𝑏}is the type of filters to use, namely Haar (2 coefficients) or Daubechies filtering with four coefficients; (2) 𝑗= {𝑅𝐺𝐵, 𝐿𝑎𝑏}is the gray level version of the RGB image or the L channel of the Lab image; and (3) 𝑘= {𝐿𝐿, 𝐴𝑙𝑙}represents if only the low–low decomposition sub-bands are considered (6 = 2 × 3 features) or all the sub-bands (24 = 2 × 3 × 4 features). We used the functions wfilters() and dwt2() of the Matlab Wavelet Toolbox. Gabor filters are sinusoidal waves modulated by a Gaussian envelope that can be used for texture classification (Randen & Husoy, 1999). The filters are applied to the images varying their frequency and orientation. Bianconi and Fernández (2007) analyzed the influence of these parameters for texture classification. After applying a set of digital Gabor filters 𝐺𝑖𝑗 (𝑥, 𝑦)with 𝑖∈ {1,…, 𝑛𝐹}and 𝑗∈ {1,…, 𝑛𝑂}, where 𝑛𝐹 and 𝑛𝑂are respectively the number of frequency and orientations, some statistical features are computed over each filtered image. We used the gabor() and imgaborfilt() functions of the Matlab Image Processing Toolbox considering the wavelengths [3,6,9,12] and orientations [0,30,60,90,120,150], recommended by Bianconi and Fernández (2007). We compute the feature vectors GaborRGB and GaborLab, both with 48 features, that include the mean and standard desviation for each filter (𝑛𝐹×𝑛𝑂= 6 × 4 = 24 filters) applied on the gray level version of the RBG image or the L channel of the Lab image, respectively. 3.3. Regression models We selected for this experimentation several state-of-art regression models of different families that provided good performances in the comparative analysis (Fernández-Delgado et al.,2019). One of them is implemented in the Octave3scientific programming language, and the remaining ones in the R4statistical computing language. Most regressors in our collection have tunable hyper-parameters, i.e., parameters whose values must be specified previously to training, that often have a strong influence on the regressor performance. In these cases, it is a good practice to try several values for each hyper-parameter in a trial-and-error procedure, and to select the value that provides the best performance on a separate data collection. This method is called ‘‘grid search’’ tuning. The following is a list of these regressors, with its tunable hyper-parameters and the list of values tried for each one in the grid search. For the regressors programmed in R, these values were provided by the getModelInfo() function of the caret R package (Kuhn,2016). 1. lm is the linear regression provided by the stats R package, which performs multivariate linear regression and has no tunable hyper-parameter (Bates & Chambers,1992). 2. svr: epsilon-support vector regression with radial basis function kernel, using the LibSVM library (Chang & Lin,2011) through its Octave interface. The regularization hyper-parameter 𝐶and 𝛾= 1∕2𝜎2, where 𝜎is the kernel spread, are tuned with values {2𝑖}15 −5 and {2𝑖}10 −10, respectively. 3. M5: regression tree (Quinlan,1992) implemented by the Weka Data Mining Software5and accessed from a R program through the RWeka package. It has no tunable hyper-parameter. 4. cubist: M5 rule-based regressor with corrections based on nearest neighbors in the training set (Quinlan,1993), implemented by the Cubist R package. Its tunable hyper-parameters are the number of neighbors [0,5,9] and the number of committees [1,10,20]. 3http://www.octave.org. 4http://www.r-project.org. 5http://www.cs.waikato.ac.nz/ml/weka. Expert Systems With Applications 206 (2022) 117765 7 E. Cernadas et al. 5. gbm: generalized boosting regression model (gbm R package) with Gaussian distribution. The tunable hyper-parameters are the maximum depth of input interactions, with values [1,2,3,4,5] and the number of trees for prediction, with values from 50 to 250 step 50. 6. rf: random forest (Breiman,2001) ensemble of averaged random regression trees (randomForest R package). The number of inputs selected at each tree (mtry) is tuned with 10 values between 2 and the number of features. 4. Experimental setup In order to achieve a fully automatic system, which can operate on-line in the meat industries, we will develop three experiments to compare the computer predictions with the sensorial procedures: 1. Experiment 1: prediction of the marbling from irregular regions of the ham muscles that are delineated and annotated by the food technology experts. In this case, the prediction is done using the same information (irregular region) as the experts. 2. Experiment 2: prediction of the marbling using square regions extracted from the irregular regions of the experiment 1. This experiment tests the reduction of performance when a region of the slice smaller than the whole irregular region is used to do the prediction. 3. Experiment 3: fully automatic extraction of square regions of interest (ROI) for each ham muscle, and prediction of the marbling using these ROIs. This experiment will test the performance loss when the ROI extraction may be suboptimal, but the method is fully automatic and does not require any food technology expert. In order to test the performance of the ROI extraction algorithm, we define the overlapping percentage (OP) as: 𝑂𝑃𝑚= 100 𝑁𝑂𝑚 𝑅𝑡 (2) where 𝑚may be biceps, semimembranosus or semitendinosus, 𝑅𝑡is the number of pixels of the ROI (in our case 𝑠2=64 ×64 =4096 pixels) and 𝑁𝑂𝑚is the number of pixels overlapped to the true region annotated by the expert for the muscle 𝑚. To test the performance of regression models in the prediction of marbling, we used the classical K-fold cross-validation methodology, which uses training and test sets. The most popular performance measures are the Pearson’s correlation coefficient (𝑅) between the true and predicted marbling, the mean absolute error (average absolute difference between the predicted and true marbling, MAE) and the root mean square error (square root of the mean squared difference between them, RMSE): 𝑅= 𝑁 ∑ 𝑖=1 (𝑦𝑖−𝑦)(𝑜𝑖−𝑜) √ √ √ √(𝑁 ∑ 𝑖=1 (𝑦𝑖−𝑦)2)(𝑁 ∑ 𝑖=1 (𝑜𝑖−𝑜)2)(3) 𝑀𝐴𝐸 =1 𝑁 𝑁 ∑ 𝑖=1 |𝑦𝑖−𝑜𝑖|(4) 𝑅𝑀𝑆𝐸 =√ √ √ √1 𝑁 𝑁 ∑ 𝑖=1 (𝑦𝑖−𝑜𝑖)2(5) where 𝑦𝑖and 𝑜𝑖are the predicted and true values of marbling respectively for ham slice 𝑖,𝑦 and 𝑜 are the mean values of {𝑦𝑖}𝑁 𝑖=1 and {𝑜𝑖}𝑁 𝑖=1, respectively, and 𝑁is the number of ham slices. The |𝑅|values can be interpreted according to Colton (1974) as: true and predicted values are not correlated at all (0–0.25), bad to moderate correlation (0.25–0.5), moderate to good (0.5–0.75), very good to excellent (0.75– 1). In our study, we use 𝐾= 4 folds or trials, devoting 𝐾− 2 = 2 Table 2 Distributions of folds in training, validation and test for each trial. 𝐾=4 Trial 1 Trial 2 Trial 3 Trial 4 Train 𝑇1={1,2} 𝑇1={2,3} 𝑇1={3,4} 𝑇1={4,1} Validation 𝑉1={3} 𝑉1={4} 𝑉1={1} 𝑉1={2} Test 𝑆1={4} 𝑆1={1} 𝑆1={2} 𝑆1={3} folds for training, one for validation and one for test. Since all the folds have the same size, 50%, 25% and 25% of the data are devoted to training, validation and test sets, respectively. In our case, each data corresponds to the image of a ham slice, and it is composed by the texture feature vector extracted from the image (input) and the marbling value (output) corresponding to that slice image. In order to guarantee that training, validation and test sets contain output values distributed across the whole range of marbling values, the data are sorted by increasing marbling. After sorting, data 1, 2, 3 and 4 are included in folds 1, 2, 3 and 4, respectively. Data 5, 6, 7 and 8 are added to folds 1, 2, 3 and 4, respectively, and so on. Let {𝑇𝑘, 𝑉𝑘, 𝑆𝑘}𝐾 𝑘=1 be the training, validation and test folds on 𝑘th trial. Trial 1 uses folds 1 and 2 for training, fold 3 for validation and fold 4 for test (see Table 2). Trial 2 uses folds 2 and 3 for training, fold 4 for validation and fold 1 for test, and analogously for trials 3 and 4. Therefore, training, validation and test sets in all the trials include data with marbling values in the whole range of values, composing training sets of higher quality that are expected to allow regressors learn better. Algorithm 2reports the experimental methodology. For each combination 𝐶𝑖of hyper-parameter values of the model, with 𝑖= 1,…, 𝑁, and for each trial 𝑘, with 𝑘= 1,…, 𝐾, the set 𝑇𝑘is used to train the model using the combination 𝐶𝑖of hyper-parameter values, while the set 𝑉𝑘is used as validation set to evaluate the performance 𝑃𝑖𝑘 of the trained model using 𝐶𝑖on 𝑉𝑘. The average 𝑃𝑖of {𝑃𝑖𝑘}𝐾 𝑘=1 is thus the performance associated to the combination 𝐶𝑖of hyper-parameter values. The process is repeated for all the combinations {𝐶𝑖}𝑁 𝑖=1, and the combination 𝐶𝐼with the highest performance 𝐼=argmax{𝑃𝑖}𝑁 𝑖=1 is selected. Then, for 𝑘= 1,…, 𝐾 the model with this best combination 𝐶𝐼is trained on the set {𝑇𝑘, 𝑉𝑘}, that includes 𝐾−1 folds, and tested on the set 𝑆𝑘(1 fold), achieving a performance 𝑃𝑘. The average of {𝑃𝑘}𝐾 𝑘=1 is the final test performance 𝑃of the model. Algorithm 2: Experimental methodology, combining 𝐾-fold cross-validation and hyper-parameter tuning. 1for 𝑖=1∶𝑁do 2for 𝑘=1∶𝐾do 3Train the model with 𝐶𝑖on dataset 𝑇𝑘 4Validate the model with 𝐶𝑖on dataset 𝑉𝑘 5𝑃𝑖𝑘=performance with 𝐶𝑖on 𝑉𝑘 6𝑃𝑖=average of {𝑃𝑖𝑘}𝐾 𝑘=1 7𝐼=argmax{𝑃𝑖}𝑁 𝑖=1 8for 𝑘=1∶𝐾do 9Train the model with 𝐶𝐼on {𝑇𝑘, 𝑉𝑘} 10 Test the model with 𝐶𝐼on 𝑆𝑘 11 𝑃𝑘=performance with 𝐶𝐼on 𝑆𝑘 12 𝑃=average of {𝑃𝑘}𝐾 𝑘=1 5. Results and discussion We present the results obtained from the different points of view: the automatic extraction of ROIs to represent each muscle (Section 5.1), the performance of marbling prediction model (Section 5.2), comparison among the performance of different prediction models (Section 5.3), the method stability with the muscles (Section 5.4) and the computational time of the different stages of the procedure (Section 5.5). Expert Systems With Applications 206 (2022) 117765 8 E. Cernadas et al. Fig. 7. Examples of square ROI extraction, in cyan the muscle contour, in pink the square ROI extracted in the centroid, in white, yellow and green the square ROI automatically extracted for BF, SM and ST muscles respectively: (a) ROI bigger than the BF muscle; (b) the automatic ROI extractor fails in SM muscle; (c) and (d) suboptimal ROIs for SM and ST, respectively. Table 3 Average percentage of overlapping pixels of the automatic extracted ROI (column Experiment 3) and the extraction from the centroid (column Experiment 2) with the true muscle annotated by the experts used in the experiment 1 for every muscle. Muscle (m) #images 𝑂𝑃𝑚 Experiment 3 Experiment 2 Biceps femoris 351 93.79 98.54 Semimembranosus 335 89.92 89.20 Semitendinosus 55 94.76 99.84 5.1. Extraction of ROIs We tested the performance of the automatic algorithm to extract square ROIs, of size 64 ×64 pixels, for each muscle described in Section 3.1, by using the measures described in Eq. (2) of Section 4 between the true muscle and the extracted ROI for each muscle 𝑚. We compared the results with the extraction done for experiment 2, in which the ROI is extracted from the centroid of the irregular region (muscle) annotated by the expert (see Section 2.3). This extraction is normally satisfactory: in experiment 2, an overlapping of 100% is achieved for the 47% of images, and the overlapping is below 80% for only the 8% of the images. In experiment 3, the extraction is perfect for the 52% of images and the overlapping is lower than 80% for only the 12% of the images. Fig. 7 shows some examples of suboptimal ROI extraction: for the ROI extraction using centroid (experiment 2), the causes of suboptimal extraction are muscles smaller than the ROI size (Fig. 7a) and muscles not rounded (Fig. 7c). For the automatic ROI extraction (experiment 3), the causes of failure are: (i) suboptimal extraction (Fig. 7d) and (ii) exchange SM and ST muscles due to artifacts in the ham slice (Fig. 7b). Table 3 shows the 𝑂𝑃𝑚achieved by the ROI extraction algorithm for each muscle 𝑚in both experiments. 5.2. Marbling prediction The texture features described in Section 3.2 are computed for all the ham slices. In the experiment 1, the texture features are computed on the irregular region, 𝑅𝑖, defined by each muscle and annotated by the food technology experts. So, Gabor and wavelet features were not computed because they must be applied on a square image. In experiment 2, squared regions 𝑅𝑆𝑖of size 64 ×64 pixels are extracted from the centroid of 𝑅𝑖and all the texture features were computed. In experiment 3, the texture features were computed on the square regions automatically extracted from the ham slice using the algorithm described in Section 3.1. Table 4 Correlation (column 𝑅) and mean absolute error (column MAE) for the marbling prediction using pure color features (upper part of the table), gray level texture features (middle part of the table) and color texture features (lower part of the table) using the support vector regression (svr) for the three experiments. Feature vector Experiment 1 Experiment 2 Experiment 3 𝑅MAE 𝑅MAE 𝑅MAE Pure color features CMRGB 0.85 0.62 0.79 0.70 0.72 0.79 CMab 0.69 0.82 0.63 0.86 0.58 0.93 CMVRGB 0.90 0.50 0.86 0.57 0.81 0.67 CMVab 0.74 0.75 0.66 0.83 0.68 0.86 FOSRGB 0.91 0.47 0.87 0.54 0.84 0.63 FOSab 0.81 0.68 0.74 0.76 0.71 0.78 Gray-level texture features HarRGB 0.91 0.49 0.88 0.56 0.80 0.72 HarLab 0.90 0.50 0.85 0.61 0.77 0.76 mlbpRGB 0.93 0.43 0.83 0.62 0.73 0.79 mlbpLab 0.92 0.45 0.81 0.66 0.71 0.80 WTrgbHaarLL – – 0.85 0.57 0.80 0.67 WTrgbHaarAll – – 0.88 0.54 0.80 0.66 WTlabHaarLL – – 0.85 0.57 0.80 0.67 WTlabHaarAll – – 0.88 0.54 0.82 0.65 WTrgbDb4LL – – 0.86 0.58 0.79 0.70 WTrgbDb4All – – 0.88 0.53 0.79 0.70 WTlabDb4LL – – 0.87 0.57 0.80 0.70 WTlabDb4All – – 0.89 0.53 0.79 0.70 GaborRGB – – 0.89 0.54 0.80 0.70 GaborLab – – 0.88 0.55 0.80 0.71 Color texture features CMVHarRGB 0.95 0.38 0.91 0.47 0.84 0.64 FOSHarRGB 0.95 0.39 0.92 0.46 0.83 0.63 CMVmlbpRGB 0.93 0.41 0.89 0.51 0.79 0.68 FOSmlbpRGB 0.94 0.40 0.90 0.48 0.80 0.67 CMVWTLabHaarAll – – 0.89 0.53 0.82 0.65 FOSWTLabHaarAll – – 0.89 0.51 0.83 0.62 CMVGaborRGB – – 0.91 0.48 0.82 0.65 FOSGaborRGB – – 0.91 0.48 0.82 0.64 Table 4 shows the correlation (𝑅) and the mean absolute error (MAE) for marbling prediction using the support vector regression (svr) for the three experiments. The feature vector FOSRGB achieved the highest 𝑅and lowest MAE among all pure color features (upper part of the table), with 𝑅=0.91, 0.87 and 0.84 in experiments 1, 2 and 3, respectively. As expected, when the performance of the automatic detection degrades, the regression accuracy decreases. But, it is still quite high for the third experiment, very good to excellent following Expert Systems With Applications 206 (2022) 117765 9 E. Cernadas et al. Fig. 8. Graphical representation of the true (blue line) and predicted (red points) ham marbling (vertical axis) for all ham slices (horizontal axis) using the svr regressor: (left panel) using the CMVHarRGB vector in the experiment 1 and (right panel) using the FOSRGB plus LBPRGB vector in the experiment 3. the Colton’s criteria. In relation with gray-level texture features (middle part of table), the use of irregular regions of the muscle (experiment 1) achieves the highest performance (𝑅=0.93 and MAE =0.43 using the texture vector mlbpRGB). In general, the performance achieved in experiment 2 is higher than the experiment 3 for all the features vectors used, noting the loss of information in the automatic selection of the ROIs. Comparing the different families of gray-level texture features for square ROIs, the wavelet features achieve the highest performance for experiment 3 (𝑅=0.82 and MAE =0.65) followed very close by the Haralick’s coefficients and Gabor features (𝑅=0.80), which means that the model is also from good to excellent. From the color space point of view, although the best results with wavelet features are achieved using the Lab color space, the difference with the use of RGB color space is not significant and the highest performance with the remaining feature vectors is better using RGB color space. We developed experiments combining the best pure color features (vectors FOSRGB and CMVRGB) with the best gray-level texture features of each texture features families. Specifically, we chose the graylevel texture vectors HarRGB for Haralick’s coefficients, mlbpRGB for local binary patterns, WTLabHaarAll for wavelet features and Gabor- RGB for Gabor features. The results are shown in the lower part of Table 4. The performance increases 0.2 for the first experiment (from 𝑅=0.93 for mlbpRGB to 𝑅=0.95 for CMVHarRGB) and 0.3 for the experiment 2 (from 𝑅=0.89 for GaborRGB to 𝑅=0.92 for FosHar- RGB), but 𝑅does not increase in experiment 3. In all the cases, the best results were achieved combining the color information (CMVRGB or FOSRGB vectors) and the gray-level texture information provided by the Haralick’s or LBP coefficients (vector HarRGB and mlbpRGB). The MAE is a performance measure easier to interpret than 𝑅from the point of view of food technology experts. Fig. 8 shows the reliability of the prediction using the svr regressor for experiments 1 (left panel) and rf regressor for experiment 3 (right panel) for the best feature vectors. The blue line represents the true marbling for each ham slice provided by the experts. The red points represent the predicted marbling by svr or rf for each ham slice. The average difference between the blue and red values for each ham slice is the MAE (0.38 and 0.60 in the left and right panels, respectively). This means that the prediction of the computer is the true value ±MAE in average. Taking in mind that the expert’s tolerance in the sensorial analysis to establish the marbling is 0.5, the prediction of the computer is comparable with the precision of the experts. 5.3. Comparing different regressors for marbling prediction In order to find the best prediction of marbling score, we applied the regressors described in Section 3.3 to the color texture features described in Section 3.2 for the three experiments. Fig. 9 shows the box plots comparing the correlation for all regressors considering all feature vectors in the experiments 1, 2 and 3 (left panel) and only in Table 5 Correlation (𝑅) for the marbling prediction using all the regressors and the best feature vector (column Feature) for the first and third experiments. Regressor Experiment 1 Experiment 3 Feature 𝑅Feature 𝑅 lm CMVmlbpRGB 0.932 ±0.003 FOSHarRGB 0.831±0.014 svr CMVHarRGB 0.948 ±0.004 CMVHarRGB 0.839 ±0.005 m5 FOSmlbpRGB 0.940 ±0.003 FOSHarRGB 0.834 ±0.017 cubist FOSmlbpRGB 0.943±0.007 FOSHarRGB 0.838±0.012 gbm FOSmlbpRGB 0.934 ±0.001 FOSRGB 0.829±0.013 rf FOSmlbpRGB 0.938±0.004 FOSmlbpRGB 0.846 ±0.008 experiment 3 (right panel). The upper and lower edges of each box indicate the 25th and 75th percentiles, respectively, the upper and lower blue segments enclose the remaining data and the red crosses are the data considered outliers. The red line inside the box is the median of the data. In both panels, there are not large difference among regressors. Considering the tree experiments (left panel), svr, rf and cubist have the highest medians, and svr has the smaller box, similar to rf but with higher median. In experiment 3 (right panel), the rf box shows the best median, followed by svr with a shorter box. Table 5 shows the color texture vector that achieved the highest correlation 𝑅, with the standard deviation over the 𝐾folds in the cross validation, for each regressor on the first (irregular regions for each muscle) and third experiments (square regions for each muscle). In experiment 3, the highest values of 𝑅are achieved by the rf regressor using the feature vector FOSmlbpRGB (𝑅= 0.846 ± 0.008). Considering that the MAE values (0.38 and 0.60) reflect the dispersion of the computer marbling prediction from the true marbling values, the computer predictions fall near the experts’ tolerance, which is 0.5. Hence, for the first experiment, the computer predicts correctly the 90.4% (67.8% for the third experiment) of samples within a tolerance of ±0.5 (the expert’s tolerance) and within a tolerance of ±1, the 99.3% (88.8% for the third experiment). 5.4. Comparing performance for different muscles In order to test if the regressor behavior is stable over the different muscles, we create two datasets with images belonging to the third experiment: (1) BFData for images of Biceps femoris (BF) muscle; and (2) SMData for images of Semimembranosus (SM) muscle. Semitendinosus muscle is not analyzed because there are few images. For BF muscle, the best performance is achieved by the feature vector FOSGaborRGB using the cubist regressor (𝑅=0.80 and MAE =0.50). For the SM muscle, the best performance is achieved by the feature vector FOSmlbpRGB using also the cubist regressor (𝑅=0.74 and MAE =0.55). Although the correlation loss is 0.12 (0.92−0.80) and 0.18 (0.92−0.74) for the BF and SM muscles respectively, the loss in MAE is only 0.04 (0.50−0.46) and 0.09 (0.55−0.46) respectively. Thus, the computer predictions with