[Question] Why are the edges of my heat map jagged and not covered?

[DeepLuckLab]

What's your question?

Thank you very much for your excellent work so far. I am currently encountering an issue: my heatmap does not cover the entire organization slice. The first image below shows the visualization of the cropped patch, and the second image shows the corresponding heatmap. I have attached the implementation code below. I would appreciate your guidance, as I am not sure where the problem lies.

  def visualize_heatmap(
          wsi,
          scores: np.ndarray,
          coords: np.ndarray,
          patch_size_level0: int,
          vis_level: Optional[int] = 2,
          cmap: str = 'coolwarm',
          normalize: bool = True,
          blur: bool = False,
          overlap: float = 0.5,
          num_top_patches_to_save: int = -1,
          output_dir: Optional[str] = "output",
  ) -> str:
      """
      Generate a heatmap visualization overlayed on a whole slide image (WSI).
  
      Args:
          wsi: Whole slide image object.
          scores (np.ndarray): Scores associated with each coordinate.
          coords (np.ndarray): Coordinates of patches at level 0.
          patch_size_level0 (int): Patch size at level 0.
          vis_level (Optional[int]): Visualization level.
          cmap (str): Colormap to use for the heatmap.
          normalize (bool): Whether to normalize the scores.
          num_top_patches_to_save (int): Number of high-score patches to save. If set to -1, do not save any. Defaults to -1.
          output_dir (Optional[str]): Directory to save heatmap and top-k patches.
  
      Returns:
          str: Path to the saved heatmap image.
      """
  
      if normalize:
          scores = rankdata(scores, 'average') / len(scores) * 100 / 100
  
      downsample = wsi.level_downsamples[vis_level]
      scale = np.array([1 / downsample, 1 / downsample])
      region_size = tuple((np.array(wsi.level_dimensions[0]) * scale).astype(int))
  
      overlay = create_overlay(scores, coords, patch_size_level0, scale, region_size)
  
      if blur:
          patch_size = np.ceil(np.array([patch_size_level0, patch_size_level0]) * scale).astype(int)
          overlay = cv2.GaussianBlur(overlay, tuple((patch_size * (1 - overlap)).astype(int) * 2 + 1), 0)
  
      img = wsi.read_region((0, 0), vis_level, wsi.level_dimensions[vis_level]).convert("RGB")
      img = img.resize(region_size, resample=Image.Resampling.BICUBIC)
      img = np.array(img)
  
      overlay_colored = apply_colormap(overlay, cmap)
      blended_img = cv2.addWeighted(img, 0.6, overlay_colored, 0.4, 0)
      blended_img = Image.fromarray(blended_img)
  
      os.makedirs(output_dir, exist_ok=True)
      heatmap_path = os.path.join(output_dir, "heatmap.png")
      blended_img.save(heatmap_path)
  
      if num_top_patches_to_save > 0:
          topk_dir = os.path.join(output_dir, "topk_patches")
          os.makedirs(topk_dir, exist_ok=True)
          topk_indices = np.argsort(scores)[-num_top_patches_to_save:]
          for idx, i in enumerate(topk_indices):
              x, y = coords[i]
              patch = wsi.read_region((x, y), 0, (patch_size_level0, patch_size_level0))
              patch.save(os.path.join(topk_dir, f"top_{idx}_score_{scores[i]:.4f}.png"))
      print('Visualize_heatmap Done!')
  
      return heatmap_path

[guillaumejaume]

@DeepLuckLab gotta check in detail, but can you try commenting out line https://github.com/mahmoodlab/TRIDENT/blob/main/trident/Visualization.py#L123C9-L123C73, thx

[DeepLuckLab]

@DeepLuckLab gotta check in detail, but can you try commenting out line https://github.com/mahmoodlab/TRIDENT/blob/main/trident/Visualization.py#L123C9-L123C73, thx

Thank you for your prompt reply. I replaced the latest Visualization.py and commented out #img = img.resize(region_size, resample=Image.Resampling.BICUBIC), but I got an error. The error message is as follows:

blended_img = cv2.addWeighted(img, 0.6, overlay_colored, 0.4, 0)
cv2.error: OpenCV(4.10.0) D:\a\opencv-python\opencv-python\opencv\modules\core\src\arithm.cpp:661: error: (-209:Sizes of input arguments do not match) The operation is neither 'array op array' (where arrays have the same size and the same number of channels), nor 'array op scalar', nor 'scalar op array' in function 'cv::arithm_op'

[DeepLuckLab]

@DeepLuckLab gotta check in detail, but can you try commenting out line https://github.com/mahmoodlab/TRIDENT/blob/main/trident/Visualization.py#L123C9-L123C73, thx

I noticed that after adding the blur code, the heatmap no longer aligns with the original image. I referenced this code from CLAM. How should I modify it? The main goal is to make the heatmap more visually appealing. Without the blur code, the heatmap can cover the tissue slice.

if blur:
patch_size = np.ceil(np.array([patch_size_level0, patch_size_level0]) * scale).astype(int)
overlay = cv2.GaussianBlur(overlay, tuple((patch_size * (1 - overlap)).astype(int) * 2 + 1), 0)

How can I obtain the blurred image while ensuring that the heatmap fully covers the original image?

[guillaumejaume]

I missed the blur argument that you took from CLAM. You confirm that, without blurring, it works fine with the original Trident code?

To make the heatmap more visually appealing, I highly recommend generating it with high overlap across patches.

eg, as done in Tuto 3:

from trident import OpenSlideWSI, visualize_heatmap
from trident.segmentation_models import segmentation_model_factory
from trident.patch_encoder_models import encoder_factory as patch_encoder_factory

# a. Load WSI to process
job_dir = './tutorial-3/heatmap_viz'
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
slide = OpenSlideWSI(slide_path='./CPTAC-CCRCC_v1/CCRCC/C3L-00418-22.svs', lazy_init=False)

# b. Run segmentation 
segmentation_model = segmentation_model_factory("hest")
geojson_contours = slide.segment_tissue(segmentation_model=segmentation_model, job_dir=job_dir, device=device)

# c. Run patch coordinate extraction
coords_path = slide.extract_tissue_coords(
    target_mag=20,
    patch_size=512,
    save_coords=job_dir,
    overlap=256, --> this is overlap in pixel from (or the stride)
)

# d. Run patch feature extraction
patch_encoder = patch_encoder_factory("conch_v15").eval().to(device)
patch_features_path = slide.extract_patch_features(
    patch_encoder=patch_encoder,
    coords_path=coords_path,
    save_features=os.path.join(job_dir, f"features_conch_v15"),
    device=device
)

#  e. Run inference 
with h5py.File(patch_features_path, 'r') as f:
    coords = f['coords'][:]
    patch_features = f['features'][:]
    coords_attrs = dict(f['coords'].attrs)

batch = {'features': torch.from_numpy(patch_features).float().to(device).unsqueeze(0)}
logits, attention = model(batch, return_raw_attention=True)

# f. generate heatmap
heatmap_save_path = visualize_heatmap(
    wsi=slide,
    scores=attention.cpu().numpy().squeeze(),  
    coords=coords,
    vis_level=1,
    patch_size_level0=coords_attrs['patch_size_level0'],
    normalize=True,
    num_top_patches_to_save=10,
    output_dir=job_dir
)

[guillaumejaume]

fixed?

[DeepLuckLab]

fixed?

I re-extracted the 256x256 patch with overlap=128. Figure 1 shows the image without blurring, and Figure 2 shows the image with blurring. However, even after blurring, the heat maps still do not overlap. I have decided not to blur the image, as this results in a normal heat map, but it may affect the aesthetics.

        def visualize_heatmap(
                wsi,
                scores: np.ndarray,
                coords: np.ndarray,
                patch_size_level0: int,
                vis_level: Optional[int] = 2,
                cmap: str = 'coolwarm',
                normalize: bool = True,
                num_top_patches_to_save: int = -1,
                output_dir: Optional[str] = "output",
                vis_mag: Optional[int] = None,
                overlay_only=False,
                blur = False,
                overlap = 0,
                filename='heatmap.png'
        ) -> str:
            """
            Generate a heatmap visualization overlayed on a whole slide image (WSI).
        
            Args:
                wsi: Whole slide image object.
                scores (np.ndarray): Scores associated with each coordinate.
                coords (np.ndarray): Coordinates of patches at level 0.
                patch_size_level0 (int): Patch size at level 0.
                vis_level (Optional[int]): Visualization level.
                cmap (str): Colormap to use for the heatmap.
                normalize (bool): Whether to normalize the scores.
                num_top_patches_to_save (int): Number of high-score patches to save. If set to -1, do not save any. Defaults to -1.
                output_dir (Optional[str]): Directory to save heatmap and top-k patches.
                vis_mag (Optional[int]): Visualization Magnification. This will overwrite vis_level
                overlay_only bool: Whenever to save the overlay only. If set to True, save the overlay on top of downscaled version of the WSI. Defaults to False.
                filename (str): file will be saved in `output_dir`/`filename`
        
            Returns:
                str: Path to the saved heatmap image.
            """
        
            if normalize:
                from scipy.stats import rankdata
                scores = rankdata(scores, 'average') / len(scores) * 100 / 100
        
            if vis_mag is None:
                downsample = wsi.level_downsamples[vis_level]
            else:
                src_mag = wsi.mag
                downsample = src_mag / vis_mag
                if not overlay_only:
                    vis_level, _ = wsi.get_best_level_and_custom_downsample(downsample)
        
            scale = np.array([1 / downsample, 1 / downsample])
            region_size = tuple((np.array(wsi.level_dimensions[0]) * scale).astype(int))
            overlay = create_overlay(scores, coords, patch_size_level0, scale, region_size)
        
            if blur:
                patch_size = np.ceil(np.array([patch_size_level0, patch_size_level0]) * scale).astype(int)
                overlay = cv2.GaussianBlur(overlay, tuple((patch_size * (1 - overlap)).astype(int) * 2 + 1), 0)
        
            overlay_colored = apply_colormap(overlay, cmap)
        
            if overlay_only:
                blended_img = overlay_colored
            else:
                img = wsi.read_region((0, 0), vis_level, wsi.level_dimensions[vis_level]).convert("RGB")
                img = img.resize(region_size, resample=Image.Resampling.BICUBIC)
                img = np.array(img)
        
                blended_img = cv2.addWeighted(img, 0.6, overlay_colored, 0.4, 0)
        
            blended_img = Image.fromarray(blended_img)
        
            os.makedirs(output_dir, exist_ok=True)
            heatmap_path = os.path.join(output_dir, filename)
            blended_img.save(heatmap_path)
        
            if num_top_patches_to_save > 0:
                topk_dir = os.path.join(output_dir, "topk_patches")
                os.makedirs(topk_dir, exist_ok=True)
                topk_indices = np.argsort(scores)[-num_top_patches_to_save:]
                for idx, i in enumerate(topk_indices):
                    x, y = coords[i]
                    patch = wsi.read_region((x, y), 0, (patch_size_level0, patch_size_level0))
                    patch.save(os.path.join(topk_dir, f"top_{idx}_score_{scores[i]:.4f}.png"))
            print('Visualize_heatmap Done')
            return heatmap_path

[guillaumejaume]

@DeepLuckLab, I'm not sure blurring is the way to go for heatmap generation. You can increase overlap to 90% of you will a real sharp heatmap.