Function extract_save_features

extract_features.py:576–642

(args)

    return np.asarray(features_), np.asarray(coords_)

def extract_save_features(args):
    # Derive the slide ID from its name.
    slide_id, _ = os.path.splitext(os.path.basename(args.input_slide))
    wip_file_path = os.path.join(args.output_dir, slide_id + "_wip.h5")
    output_file_path = os.path.join(args.output_dir, slide_id + "_features.h5")

    os.makedirs(args.output_dir, exist_ok=True)

    # Check if the _features output file already exist. If so, we terminate to avoid
    # overwriting it by accident. This also simplifies resuming bulk batch jobs.
    if os.path.exists(output_file_path):
        raise Exception(f"{output_file_path} already exists")

    # Open the slide for reading.
    wsi = openslide.open_slide(args.input_slide)

    # Decide on which slide level we want to base the segmentation.
    seg_level = wsi.get_best_level_for_downsample(64)

    # Run the segmentation and  tiling procedure.
    start_time = time.time()
    tissue_mask_scaled = create_tissue_mask(wsi, seg_level, method=args.method)
    filtered_tiles = create_tissue_tiles(wsi, tissue_mask_scaled, args.tile_size)

    # Build a figure for quality control purposes, to check if the tiles are where we expect them.
    qc_img = make_tile_QC_fig(filtered_tiles, wsi, seg_level, 2)
    qc_img_target_width = 1920
    qc_img = qc_img.resize((qc_img_target_width, int(qc_img.height / (qc_img.width / qc_img_target_width))))
    qc_img_file_path = os.path.join(args.output_dir, f"{slide_id}_features_QC.png")
    qc_img.save(qc_img_file_path)
    print(f"Finished creating {len(filtered_tiles)} tissue tiles in {time.time() - start_time}s")

    # Save QC figure.
    qc_img_file_path = os.path.join(
        args.output_dir, f"{slide_id}_N{len(mergedpatches)}mergedpatches_distThreshold{args.dist_threshold}_corrThreshold{args.corr_threshold}.png"
    )

    # Extract the rectangles, and compute the feature vectors. Example using EsVIT. 
    device = torch.device("cuda") 
    model, _, depths = load_encoder_esVIT(args, device)
    
    features, coords = extract_features(
        model,
        device,
        wsi,
        filtered_tiles,
        args.workers,
        args.out_size,
        args.batch_size,
        n_last_blocks = args.n_last_blocks, 
        avgpool_patchtokens = args.avgpool_patchtokens,
        depths = depths,
    )
    
    print(f'Number of features N={len(features)}')
    # Merging nearby patches with similar semantic. 
    mergedpatches = mergedpatch_gen(features, coords, dist_threshold=args.dist_threshold, corr_threshold=args.corr_threshold)
    print(f'Merging step => N={len(mergedpatches)}')