3.1.1.

Sinusoid extraction

We first apply an orientation-selective (OS) filter to improve the contrast of the sinusoid boundaries. We then apply the expectation–maximization (EM) algorithm to cluster bright pixels. Figure 4(a) shows the original image, and Fig. 4(b) shows the result of sinusoid extraction. The effectiveness of the OS filter is shown in Ishikawa et al.,¹⁵ where details of this step can be found.

Fig. 4

(a) Original HE-stained liver tissue image ($1024\times 1024$) and (b) result of sinusoid extraction. In (b), the white pixels indicate the sinusoid regions and the black pixels indicate the other regions.

3.1.2.

Stroma extraction

Stroma extraction consists of two steps. In the first step, we calculate the fiber probability of each pixel using color and texture features. The color and texture features are calculated over the surrounding $64\times 64$ region of each pixel. Color features consist of the (1) average, (2) variance, (3) skewness, and (4) kurtosis acquired from the histogram of each color channel, which corresponds to $4\text{features}\times 3\text{color channels}=12\text{features}$. The texture features consist of the (1) angular second moment, (2) contrast, (3) correlation, and (4) homogeneity of Hararick’s texture features, calculated from the gray-level co-occurrence matrix (GLCM). The pixel brightness values are normalized to the range of 0 to 1 and are quantized into eight levels.

For the GLCM-based features, extraction was performed using a four-direction GLCM (0, 45, 90, and 135 deg) with a distance of 1. The four statistical quantities, calculated from the four co-occurrence matrices using different direction offsets, are summarized by two quantities, the mean and variance. Thus, the final texture feature is a 24-dimensional feature ($\text{mean}\times \text{variance}\times 4\text{statistical features}\times 3\text{color channel}=24\text{features}$). The fiber probability is calculated for each pixel using a linear support vector machine based on the 36-dimensional color and texture features. We constructed our learning data set using 4000 image patches of $64\times 64\text{pixels}$, half of which were fiber regions. The remaining half included no fiber at all. An example of pixel-wise fiber probability indicated by pseudocolor is shown in Fig. 5(a). Figure 5(b) shows the stroma area in green, extracted by a threshold of 0.3. Basically, we have successfully extracted the stroma area. However, overdetection has occurred in the cytoplasm and underdetection has occurred in the stroma. A larger amount of stroma has been detected in the border between the stroma and the cytoplasm, and hence the border is not accurate. To solve this problem, we apply the next step.

Fig. 5

(a) Fiber probability, (b) extraction result using fiber probability ($\text{threshold}=0.3$), (c) SLIC result, (d) number of lymphocytes per superpixel, (e) mean fiber probability per superpixel, and (f) automatically extracted stroma region (indicated by green). For images (a), (d), and (e), higher values are represented by warm colors, and lower probabilities are represented by cool shades.

In the second step, we extract stroma based on the “superpixel” image processing technique. We generated superpixels using simple linear iterative clustering (SLIC). SLIC is a fast and memory-efficient adaptation of the K-means algorithm to generate superpixels. Details can be found in Ref. 20. Figure 5(c) shows the SLIC result. Stroma identification is performed per superpixel using fiber probability and the number of lymphocytes. We empirically determined the thresholds for fiber probability and number of lymphocytes. Specifically, we determined the optimal threshold by applying various thresholds to a large number of images. The purpose of stroma extraction is to improve the trabeculae extraction. Overdetection in stroma extraction is a bigger problem than underdetection. Therefore, superpixels with a mean fiber probability of more than 0.3 are extracted as stroma. Additionally, superpixels containing more than seven lymphocytes are extracted as stroma. Lymphocytes can be extracted by template matching. In this paper, we prepared a circle template with a diameter of $7.3\mu \mathrm{m}$. Figure 5(d) shows a pseudocolor representation of the number of lymphocytes per superpixel, and Fig. 5(e) shows the pseudocolor representation of mean fiber probability per superpixel. Higher probability is indicated by warm colors, and lower probability is indicated by cool shades. Figure 5(f) shows the extracted stroma area in green. Compared with Fig. 5(b), the errors of Fig. 5(f) with respect to overdetection and the boundaries of the stroma and cytoplasm are reduced.

3.2.1.

Preprocessing

The lower-left boundary of the large fat droplet shown in the center of Fig. 7(a) has interrupted parts on its boundary. Such an interrupted boundary, which occurs on some droplets, affects the results of the circularity calculation. OS filtering is a sinusoid extraction preprocessing step because it connects the boundaries of a sinusoid. For the same reason, the OS filter is effective for fat droplet extraction. Therefore, we adopt an OS filter to smooth the fat droplet boundaries (step 1). The result of the OS filter is shown in Fig. 7(b). To clarify its effects, Figs. 7(a) and 7(b) were binarized at a threshold of 175 and the results are shown in Figs. 7(c) and 7(d), respectively. It can be observed that the boundaries are clear when the OS filter is applied.

Fig. 7

Results of OS filter: (a) fat droplets, (b) result of OS filter, (c) binary image of (a), and (d) binary image of (b).

3.2.2.

Extracting candidate fat droplet regions

We next extract the candidate fat droplet regions using an EM algorithm (step 2). The fat droplet basically looks whiter than the sinusoid because its fat and lipids liquidate during the process of creating a paraffin section for the slide, leaving a completely hollow space inside. Therefore, we classify the RGB image into 10 classes using the EM algorithm and extract images in the brightest class as fat droplet candidate regions.

3.2.3.

Specific of fat droplets

We next label the results of extraction and distinguish fat droplets from sinusoids and other structure in each labeled area. Fat candidate areas are classified by a random forest (step 3). As to feature values, circularity, eccentricity, area, complexity, clockwise turning ratio, outline concentration ratio, turning angle median value, derivative matching feature, fractal dimension, number of edge pixels, contrast, intensity mean, intensity variance, outside intensity mean, outside variance of intensity, outside contrast, number of nuclei, area change ratio, and outside number of nuclei were used.

Figure 8(a) shows the final detected fat droplet regions. Finally, we modify trabeculae extraction results using the fat extraction results (step 4). We correct the cytoplasm results by excluding the fat droplet areas extracted from the sinusoidal areas. To be more specific, the fat areas of Fig. 8(a) are used to exclude the detected fat droplets in the sinusoidal area of Fig. 8(b). The corrected results using this proposed method is shown in Fig. 8(c).

Fig. 8

Sinusoid extraction in tissue that includes fat droplets: (a) result of fat droplet extraction, (b) result of trabecular extraction after extracting the fibers and sinusoids, and (c) result of trabecular extraction after extracting the fibers, sinusoids, and fat droplets. In (a)–(c), the white pixels indicate the target regions (fat droplets and trabeculae) and the black pixels indicate the other regions.

5.1.1.

Nuclear–cytoplasmic ratio for each cord label

Figure 16(a) shows a boxplot of the “mean nucleus-to-cytoplasm ratio” of Huang et al.,⁹ and Fig. 16(b) shows a boxplot of the N/C ratio distribution calculated using the method explained in Sec. 4. To calculate the result shown in Fig. 16(b), the N/C ratio of the per-unit area ($2174\times 2174\text{pixels}$) was derived. The median N/C ratio for each cord label ($p$-value $<0.01$) and mean nucleus-to-cytoplasm ratio ($p$-value $<0.01$) were evaluated by the Kruskal–Wallis test. Figures 16(a) and 16(b) present the result of a multiple comparison test, which shows that the difference between the cancer and noncancer labels was significant. To verify the separation of BG and G1 by the proposed method, we computed the receiver operating characteristic (ROC) curve for classification by a logistic regression of the N/C ratio and proposed method [see Fig. 16(c)]. The area under the curve (AUC) was significantly higher for the median N/C ratio for each cord label than it was for the mean nucleus-to-cytoplasm ratio ($p$-value $<0.05$). The median of the N/C ratio calculated for all cords is more clearly separated than the mean nucleus-to-cytoplasm ratio calculated using the total areas of nuclear and cytoplasmic regions. This indicates that the proposed method improves the separation of well-differentiated HCC from noncancer tissue. Table 1 shows the AUC of the N/C ratio and median of N/C ratio for each cord label.

Fig. 16

Boxplot of the (a) N/C ratio (Huang et al.), (b) the median of N/C ratio of each cord label, and (c) ROC curve.

Table 1

AUC and p-value for each feature.

5.1.2.

Nuclear density for each cord label

Figure 17(a) shows a boxplot of the median nuclear density method of Kiyuna et al.,¹³ and Fig. 17(b) shows a boxplot of the nuclear density distribution calculated using the method explained in Sec. 4. The median nuclear density of each cord label ($p$-value $<0.01$) and nuclear density of Kiyuna et al. ($p$-value $<0.01$) were evaluated by the Kruskal–Wallis test. Figures 17(a) and 17(b) present the result of a multiple comparison test, which shows that the difference between the cancer and noncancer labels was significant. To verify the separation of BG and G1 by the proposed method, we computed the ROC curve for classification using a logistic regression of the nuclear density and proposed methods [see Fig. 17(c)]. An ROC curve was obtained for the median nuclear density and median nuclear density distribution for each cord label threshold. In this figure, the true positive rate is plotted on the $y$-axis versus the false positive rate on the $x$-axis. There was no significant difference between nuclear density of the proposed method and Kiyuna et al.’s method ($P=0.836$). Table 1 shows the AUCs of the nuclear density and median of nuclear density of each cord label. The method of Kiyuna et al. calculates nuclear density per unit area. The proposed method calculates nuclear density per unit trabeculae. Therefore, if the image includes substances other than cells, like stroma or regions of glass, the results of the proposed and conventional methods would be different. No significant difference was observed in this experiment because an image whose entire space was occupied by cells was used.

Fig. 17

Boxplot of the (a) nuclear density, (b) median of nuclear density for each cord label, and (c) ROC curve.

5.1.3.

Number of layers for each cord label

Figure 18(a) shows a boxplot of the trabecular thickness of Kiyuna et al.,¹³ and Fig. 18(b) shows a boxplot of the distribution of the number of layers, calculated using the method explained in Sec. 4. The median number of layers of each cord label ($p$-value $<0.01$) and median trabecular thickness feature of Kiyuna et al. ($p\text{-}\text{value}=0.478$) were evaluated by the Kruskal–Wallis test. Figures 18(a) and 18(b) present the result of a multiple comparison test, which shows that the difference between the cancer and noncancer labels was significant. To verify the separation of BG and G1 by the proposed method, we computed the ROC curve for classification by a logistic regression of the trabecular thickness and proposed method [see Fig. 18(c)]. An ROC curve was obtained for trabecular thickness and median number of layer distributions for each cord label threshold. In this figure, the true positive rate is plotted on the $y$-axis versus the false positive rate on the $x$-axis. There was no significant difference between the nuclear density of the proposed method and Kiyuna et al.’s method ($P=0.141$). Table 1 shows the AUCs for trabecular thickness and the median of the number of layers of each cord label. However, the number of layers of the Kruskal–Wallis test had a significant difference for BG and G1 as well as for BG and G2. Therefore, to verify the separation of BG and G1,G2 by the proposed method, we computed the ROC curve for classification by a logistic regression of the trabecular thickness and proposed method (see Fig. 19). The AUC was significantly higher for the median N/C ratio of each cord label than it was for the mean nucleus-to-cytoplasm ratio ($p$-value $<0.01$).

Fig. 18

Boxplot of the (a) trabecular thickness, (b) number of layers for each cord label, and (c) ROC curve.

Fig. 19

ROC curve of trabecula thickness and number of layers of BG and G1,G2.

The trabecular thickness method of Kiyuna et al. is not highly accurate. Although the trabecular thickness is a classification feature for HCC, it is difficult to diagnose HCC using trabecular thickness alone. Therefore, Kiyuna et al. proposed the use of a combination of various features such as nuclear density and NC ratio. The proposed method obtains a higher accuracy than the trabecular thickness method. It is believed that classification based on the number of layers alone is difficult because the structural changes in G1 compared to normal tissue are very small. As the proposed method is an approach that quantifies the number of layers, which is just one part of the diverse structural changes that occur in liver cells, it is not useful for all situations. Therefore, it is believed that the number of layers would be best used in combination with various other features such as nuclear density and NC ratio. In fact, in a study conducted by Aziz et al.,²³ the accuracy of the difficult classification of G1 was improved from 65% to 77% by combining nuclear features with the proposed method.