Deep Convolutional Neural Networks for Breast Cancer Histology Image Analysis

In plain words

Pathologists look at thin tissue slides, stained pink and purple, to tell whether a breast biopsy is normal, harmless (benign), confined to the milk duct (in situ) or spreading into the surrounding tissue (invasive). Their readings often disagree, so computer help is useful. This paper builds such a helper for four-way sorting of the slides.

Training a large image network from scratch needs far more pictures than the 400 available here. The authors therefore use networks already trained on ordinary photos as fixed measuring tools. Each slide is recoloured and cropped many times, every crop becomes a list of numbers, and those numbers are summarised per slide. A simple but strong classifier, gradient-boosted trees, learns from the summaries. Several networks and crop sizes are combined by averaging their predictions.

The lesson the authors draw is that many varied views of a small dataset, plus averaging over models, is what makes a few hundred images workable here.

TL;DR

• Frozen ImageNet encoders (ResNet-50, InceptionV3, VGG-16) turn heavily augmented H&E crops into one pooled vector per augmentation.
• LightGBM classifiers are trained per encoder and crop size, with 10 stratified folds and 5 seeds, and their probabilities are averaged.
• Author-reported: 87.2% four-class accuracy and 93.8% carcinoma-vs-non-carcinoma accuracy under 10-fold cross-validation on the ICIAR 2018 data.
• Reproduction needs the committed fold file and the code's LightGBM values, because the paper leaves those out.

Problem

Breast biopsy classification by pathologists is tedious and inconsistent between readers §1. The training set is small: 400 images across four classes §2.1. The authors report that fine-tuning a pretrained network did not work well for them, but give no numbers for that comparison §2.2. Their approach assumes that an ImageNet-pretrained network's features plus strong augmentation can stand in for task-specific training §2.2.

Method

The pipeline has five stages. Stages 1 to 4 run once per image and encoder; stage 5 is the classifier.

1. Data. 400 H&E images of 2048×1536 px at 200× magnification, 0.42 µm per pixel, in four balanced classes: normal, benign, in situ carcinoma, invasive carcinoma §2.1 fig. 1. The labels come from two medical experts §2.1.
2. Stain normalisation. Normalise the amount of H&E stain on each image §2.3. The paper cites Macenko et al. for this. The code instead implements a variant from Araújo et al. with fixed constants feature_extractor.py.
3. Downscale. Halve the image to 1024×768 px §2.3. The code uses only this scale feature_extractor.py.
4. Colour augmentation. Make 50 augmentations per image. Each one decomposes the RGB colour into haematoxylin and eosin (H&E) components (Ruifrok and Johnston), then multiplies the H and E magnitude of every pixel by two independent uniform draws from [0.7, 1.3] §2.3. The code fixes the decomposition matrices and leaves the third (residual) channel unscaled feature_extractor.py.
5. Crops. From each augmented image take 20 random square crops, each 400 or 650 px §2.3. The code also rotates each crop by a random multiple of 90° and applies random flips feature_extractor.py. The paper does not mention these.
6. Encoders. Use ImageNet-pretrained ResNet-50 and InceptionV3 with global average pooling, giving 2048 features each §2.4. For VGG-16, average-pool the outputs of four convolutional layers from blocks 2 to 5 (128, 256, 512 and 512 channels), concatenated to 1408 features §2.4 fig. 3. The code names the layers as block2_conv2, block3_conv3, block4_conv3 and block5_conv3 models.py. Neither the paper nor the code shows a final normalisation step before the classifier.
7. Pooling. Pool the 20 crop vectors of one augmentation into a single vector with a 3-norm, applied element by element §2.3 eq. 1. The code uses the same power mean with p = 3 feature_extractor.py.
8. Descriptors. Each image therefore has 50 descriptors per encoder and crop size. Across 3 encoders and 2 crop sizes that is 300 §2.3.
9. Classifier. For each encoder and crop size, train one multiclass LightGBM model per fold and seed. Each training row is one augmented descriptor, so an image contributes 50 rows §2.5 train_lgbm.py.
10. Folds and seeds. Split at the image level into 10 stratified folds, so all 50 descriptors of an image stay in one fold §2.5 utils.py. Train 5 seeded models per fold §2.5 train_lgbm.py.
11. Blend. Average the predicted probabilities over augmentations, seeds and the six encoder/crop models, then take the argmax §2.5 crossvalidate_blending.py. Test predictions use all fold models for each encoder and crop size submission.py.
12. Two-class task. The paper reports a carcinoma-vs-non-carcinoma task but does not say how it is derived §3. The code takes the four-class blend, maps argmax labels to two classes, and builds two-class scores by taking the maximum over {normal, benign} and over {in situ, invasive} crossvalidate_blending.py. It does not train a separate binary model.

Pooling equation, as written in the paper eq. 1:

\mathbf{d}_{pool}=\left(\frac{1}{N}\sum\limits_{i=1}^N(\mathbf{d}_i)^p\right)^{\frac{1}{p}} eq. 1

with p = 3, N = 20 crops per augmentation, \mathbf{d}_i the crop descriptor and \mathbf{d}_{pool} the pooled descriptor §2.3.

Pseudocode for stages 3 to 8. The paper does not give an algorithm box, so this one is ours:

# Descriptors for one image  [§2.3] [§2.4]
# Input:  image (H, W, 3) RGB, one encoder, crop size s in {400, 650}
# Output: descriptors (50, D), one row per colour augmentation

def image_descriptors(image, encoder, crop_size):
    # Downscale by 0.5, then normalise the stain once  [§2.3]
    image = normalize_stain(downscale(image, 0.5))

    descriptors = []
    for augmentation in range(50):
        # Scale H and E by two random factors in U(0.7, 1.3)  [§2.3]
        augmented = scale_hematoxylin_eosin(image, 0.7, 1.3)

        # 20 random square crops of size s per augmentation  [§2.3]
        crops = random_crops(augmented, crop_size, 20)

        # Encode the crops, then pool them with a 3-norm  [eq. 1]
        features = encoder(crops)                       # (20, D)
        descriptors.append(power_mean(features, 3))     # (D,)

    return stack(descriptors)                           # (50, D)

# Cross-validated blend  [§2.5]
# Input:  descriptors per model, image labels, 10 folds
# Output: predicted class per image

def cross_validated_blend(descriptors, labels, folds):
    # A held-out image is scored only by models that never saw it  [§2.5]
    scores = {}
    for fold_index, fold in enumerate(folds):
        for model in models:                     # encoder x crop size
            for seed in range(5):
                # Train on the other nine folds  [code: train_lgbm.py#L87]
                booster = train_boosted_trees(
                    rows=rows_of(fold.train, model),    # 50 per image
                    labels=labels_of(fold.train),
                    rounds=70,
                    seed=1 + seed * 10 + fold_index,
                )
                # Keep 50 probabilities per held-out image
                scores[model, fold_index, seed] = booster.predict(
                    rows_of(fold.test, model)
                )

    # Average over seeds, augmentations and models  [§2.5]
    blended = mean_over_seeds_augmentations_models(scores)
    return argmax(blended)

Hyperparameters

The paper names LightGBM and the encoders but gives almost no values §2.4 §2.5. Values below come from the code, which means they are not stated in the paper.

Results

Author-reported. Setup: ICIAR 2018 Grand Challenge data, 400 H&E images, four balanced classes, 10-fold stratified cross-validation §2.1 §3. The paper gives no numeric comparison with other methods. The claim that the approach beats prior automated analyses §4 is not backed by a table.

Full per-fold table as printed in the paper tab. 1. Accuracy in percent; "Model fusion" is the blend described in stage 11.

tab. 1

Headline numbers:

The confusion matrix is in fig. 4b; the paper does not say which setpoint it uses. The stated 9 + 5 = 14 misses give 186/200 = 93.0%, which matches the 0.50 sensitivity. This is our arithmetic, not a stated fact.

Paper vs code

Code at commit 3eedb3c2de946e30ab207389a649d2a018109f7d, owned by the first author, MIT LICENSE.

• Number of models. Paper says 600 models, from 10 folds × 5 seeds × 4 scale-and-crop combinations × 3 encoders §2.5. Code runs only scale 0.5, so it trains 10 × 5 × 2 × 3 = 300 train_lgbm.py. The paper's pipeline also uses only the downscaled images §2.3.
• Crop augmentation. Paper describes random crops §2.3. Code rotates each crop by a random multiple of 90° and flips it feature_extractor.py. A free-rotation variant is written but commented out feature_extractor.py.
• Stain normalisation. Paper cites Macenko et al. §2.3. Code implements the Araújo et al. variant with fixed constants feature_extractor.py.
• Colour augmentation. Matches the paper: two independent uniform factors on H and E, and the residual channel is unchanged §2.3 feature_extractor.py.
• VGG-16 layers. Paper says blocks 2 to 5 §2.4. Code names one convolutional layer per block models.py. Channel counts match.
• Augmentations and crops. Paper says each image is "represented by 20 crops" and then describes 50 augmentations §2.3. Code takes 20 crops for each of the 50 augmentations feature_extractor.py.
• Test-time descriptors. Paper says 300 descriptors per test image §2.5. README says 50 per encoder and crop size [README.rst#L61]. Code produces 50 per encoder and crop size and averages over all 6, which gives 300 in total submission.py.
• Blending. Paper says "averaged over all augmentations and models" §2.5. Code does that crossvalidate_blending.py. A logistic-regression blender also exists crossvalidate_blending.py but the reported numbers use the plain average.
• Two-class results. Paper reports them without describing the derivation §3. Code derives them from the four-class blend, using the maximum over each pair of classes crossvalidate_blending.py. The commented-out alternative sums instead crossvalidate_blending.py. Code dumps ROC scores but computes neither AUC nor the sensitivity and specificity at the two setpoints crossvalidate_blending.py. Those values cannot be reproduced from the repo as shipped.
• Recorded outputs. Code docstring records a four-class blend with mean 0.873 and std 0.0261, and a two-class blend with mean 0.938 and std 0.023 crossvalidate_blending.py. The paper prints 87.2 and 93.8. The fold values average to 87.25, so the paper rounds down.
• LightGBM values. Not in the paper §2.5. Code values are in the hyperparameter table. The bagging_fraction in the log comments (0.66) differs from the active parameter (0.69) train_lgbm.py.

Pitfalls

• LightGBM values are missing from the paper. Use the code's values, as listed in the hyperparameter table. The code also has no early stopping: the test fold is passed as a validation set, but the number of rounds is fixed at 70 train_lgbm.py.
• Bagging is off. bagging_freq is 0, so bagging_fraction has no effect in LightGBM train_lgbm.py.
• Folds cannot be regenerated exactly. The committed data/folds-10.pkl is needed to match the reported folds. The fold generator shuffles with no fixed seed utils.py. Regenerating folds gives different results.
• Augmentation randomness is not seeded. hematoxylin_eosin_aug builds RandomState(seed) with seed = None feature_extractor.py, so it ignores the global np.random.seed(13) set at import feature_extractor.py. Crop sampling does use the global RNG feature_extractor.py. Exact descriptors will not reproduce across runs.
• "×300" is not the size of one model's training set. The paper's ×300 counts augmentations across all encoder and crop combinations §2.5. Each LightGBM model sees 50 rows per training image, from one encoder and one crop size train_lgbm.py.
• Fusion gain is smaller than the text says. The paper says the fused model is 4 to 5 points better than any constituent §3. Table 1 gives 87.2 against individual means from 82.5 to 84.2, which is a gap of 3.0 to 4.7 points tab. 1. The paper also says the blend's standard deviation is about half the models' average, and the table gives 2.6 against a mean of about 4.8 tab. 1.
• Fine-tuning comparison is unsupported. The paper says fine-tuning did not work, with no data shown §2.2.
• No baseline table. Claims of superiority to earlier work §1 §4 are not backed by a numeric comparison in this paper.
• Data is not in the repo. Obtain the ICIAR 2018 images under the challenge terms; the README says the data cannot be redistributed [README.rst#L131]. Images must be .tif files feature_extractor.py.
• Software is pinned to 2018. The README lists Keras and Theano, and says TensorFlow was not tested [README.rst#L95]. The repo has no requirements file, so versions are not pinned.
• Test features need explicit paths. Feature extraction defaults to data/train, so test features need the --images and --features flags feature_extractor.py.

Lineage

• Builds on: ImageNet-pretrained ResNet-50, InceptionV3 and VGG-16 (Keras) §2.4; LightGBM §2.2; deep convolutional feature pooling §2.3.
• Compared against: the paper names prior breast histology work (Spanhol et al. 2016; Araújo et al. 2017; Bejnordi et al. 2017) and says it outperforms them, without a numeric comparison §4.
• Superseded by: not known from wiki content.