nuclei_process.py

#!/usr/bin/env python

# -*- coding: utf-8 -*-
"""GA_3rd.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1lWkTNoLKKu5xCC4QbzamJH6eI-VC3l-U
"""

import json
import numpy as np
import matplotlib.path as mplPath
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import scipy.ndimage
import os
from shapely.geometry import box

# hyper parameters
chrom_train_path = '/work/home/acvwd4uw3y181/rsliu/data000/Chromosome/train'
chrom_train_json_path = '/work/home/acvwd4uw3y181/rsliu/data000/Chromosome/train.json'
isbi_train_path = '/work/home/acvwd4uw3y181/rsliu/data000/ISBI_detection/isbi_train'
isbi_train_json_path = '/work/home/acvwd4uw3y181/rsliu/data000/ISBI_detection/isbi_train_v3.json'
kaggle_train_path = '/work/home/acvwd4uw3y181/rsliu/data000/kaggle2018/stage_1_train'
kaggle_train_json_path = '/work/home/acvwd4uw3y181/rsliu/data000/kaggle2018/train.json'
CN_train_path = '/work/home/acvwd4uw3y181/rsliu/data000/cluster_nuclei/train'
CN_train_json_path = '/work/home/acvwd4uw3y181/rsliu/data000/cluster_nuclei/train.json'

# hyper parameters
training_json_path = CN_train_json_path
training_image_path = CN_train_path + '/'

window_size = 8
GA_factor = 0.5
overlap_HL = 0.1

if not os.path.exists(training_image_path + f'anomap_ws{window_size}_f{GA_factor}_ol{overlap_HL}'):
    os.makedirs(training_image_path + f'anomap_ws{window_size}_f{GA_factor}_ol{overlap_HL}', exist_ok=True)
weighted_image_path = training_image_path + f'anomap_ws{window_size}_f{GA_factor}_ol{overlap_HL}'

target_category_id = 1

def show_weighted_image(weighted_image, original_image):
  fig, axs = plt.subplots(1, 3, figsize=(15, 5))
  axs[0].imshow(original_image, cmap='gray')
  axs[0].set_title('Original Image')
  axs[1].imshow(weighted_image, cmap='hot')
  axs[1].set_title('Weighted Image')
  axs[2].imshow(original_image, cmap='gray')
  axs[2].imshow(weighted_image, cmap='hot', alpha=0.5)  # alpha参数设置透明度
  axs[2].set_title('Overlay Image')
  plt.show()

def save_weighted_image(weighted_image, original_image, filename, save_dir):
    fig, axs = plt.subplots(1, 3, figsize=(15, 5))
    axs[0].imshow(original_image, cmap='gray')
    axs[0].set_title('Original Image')
    axs[1].imshow(weighted_image, cmap='hot')
    axs[1].set_title('Weighted Image')
    axs[2].imshow(original_image, cmap='gray')
    axs[2].imshow(weighted_image, cmap='hot', alpha=0.5)  # alpha参数设置透明度
    axs[2].set_title('Overlay Image')

    # 保存图像到文件
    plt.savefig(os.path.join(save_dir, filename + '_weighted.png'))
    plt.close(fig)

def save_weighted_npy(weighted_image, filename, save_dir):
    # 保存NumPy数组到.npy文件
    np.save(os.path.join(save_dir, filename + '_weighted.npy'), weighted_image)

def show_weighted_images_per_instance(weighted_images, original_image):
    num_instances = len(weighted_images)
    #print(num_instances)
    fig, axs = plt.subplots(num_instances, 3, figsize=(15, 5*num_instances))
    if num_instances == 1:
        axs = np.expand_dims(axs, axis=0)
    for i in range(1):
        axs[i, 0].imshow(original_image, cmap='gray')
        axs[i, 0].set_title('Original Image')
        axs[i, 1].imshow(weighted_images[i], cmap='hot')
        axs[i, 1].set_title('Weighted Image for Instance {}'.format(i+1))
        axs[i, 2].imshow(original_image, cmap='gray')
        axs[i, 2].imshow(weighted_images[i], cmap='hot', alpha=0.2)  # alpha参数设置透明度
        axs[i, 2].set_title('Overlay Image for Instance {}'.format(i+1))

    plt.tight_layout()
    plt.show()

def save_weighted_npy_per_instance(weighted_images, save_dir, instances_ids):
    for i, weighted_image in enumerate(weighted_images):
        # 使用索引作为实例的id，生成文件名
        filename_instance = "{}.npy".format(instances_ids[i])
        # 保存NumPy数组到.npy文件
        np.save(os.path.join(save_dir, filename_instance), weighted_image)

def find_masks(num_instances, this_image, instances_in_this_image, image_shape):
  masks_in_this_image = []
  rows, cols = image_shape[0], image_shape[1]
  for i in range(num_instances):
    if len(instances_in_this_image[i]["segmentation"]) == 1:
      mask = np.zeros(image_shape[0:2])
      mask_coordinates = np.array(instances_in_this_image[i]['segmentation']).reshape(-1, 2)
      mask_path = mplPath.Path(mask_coordinates)
      y, x = np.mgrid[:rows, :cols]
      mask[mask_path.contains_points(np.vstack((x.flatten(), y.flatten())).T).reshape(rows, cols)] = 1
      masks_in_this_image.append(mask)
    else:
      masks = []
      for seg in instances_in_this_image[i]["segmentation"]:
          mask = np.zeros(image_shape[0:2])
          # 将序列reshape为(-1,2)
          mask_coordinates = np.array(seg).reshape(-1,2)
          # 创建路径
          mask_path = mplPath.Path(mask_coordinates)
          # 创建网格
          y, x = np.mgrid[:rows, :cols]
          # 更新掩码
          mask[mask_path.contains_points(np.vstack((x.flatten(), y.flatten())).T).reshape(rows, cols)] = 1
          #print(np.sum(mask>0))
          masks.append(mask)
      mask = np.zeros(image_shape[0:2])
      for m in masks:
        mask += m
      mask = mask>0
      masks_in_this_image.append(mask)
  return masks_in_this_image

def find_grad_fields(instances_masks, factor):
  grad_fields = {"grad_magnitude":[], "grad_direction":[]}
  for mask in instances_masks:
    rows, cols = mask.shape[0], mask.shape[1]
    y1, x1 = np.where(mask == 1)[0],np.where(mask == 1)[1]

    center_x, center_y = np.mean(x1), np.mean(y1)

    z = np.zeros_like(mask, dtype = float)
    for m in range(rows):
      for n in range(cols):
        z[m, n] = (m - center_y)**2 + (n - center_x)**2

    grad_y, grad_x = np.gradient(z)
    gradient_magnitude = np.sqrt(grad_y**2 + grad_x**2)
    gradient_direction = np.arctan2(grad_y, factor * grad_x)
    gradient_direction[mask==0] = -4
    grad_fields["grad_magnitude"].append(gradient_magnitude)
    grad_fields["grad_direction"].append(gradient_direction)

  return grad_fields

def find_GA_per_instance_original(image_info, window_size, sliding_range):
  weighted_images_original = []
  rows, cols = image_info["image_shape"][0],image_info["image_shape"][1]
  for ins in range(1):
    if image_info["instances_categories"][ins] != target_category_id:
      weighted_image = np.zeros((rows, cols))
    else:
      diff_array = np.zeros((rows-window_size+1, cols-window_size+1))
      # print('slide y from to:',max(0, sliding_range[1]), min(rows-window_size+1, sliding_range[1]+sliding_range[3]))
      # print('slide x from to:',max(0, sliding_range[0]), min(cols-window_size+1, sliding_range[0]+sliding_range[2]))
      for i in range(max(0, sliding_range[1]), min(rows-window_size+1, sliding_range[1]+sliding_range[3])):
        for j in range(max(0, sliding_range[0]), min(cols-window_size+1, sliding_range[0]+sliding_range[2])):
          mask_window = image_info["masks"][ins][i:i+window_size, j:j+window_size]
          grad_window = image_info["grad_fields"]["grad_direction"][ins][i:i+window_size, j:j+window_size]
          if np.any(mask_window):
            diffs = []
            for k in range(image_info["num_instances"]):
              if k != ins and image_info["instances_categories"][k]==target_category_id:
                diff1 = np.std(grad_window)
                window_1D_k = grad_window.reshape(-1,)
                window_1D_k = np.delete(window_1D_k, np.where(window_1D_k==-4))
                diff2 = np.std(window_1D_k)
                diffs.append(min(diff1, diff2))
            diff_array[i,j] = np.mean(diffs)
      weighted_image = np.zeros((rows,cols))
      for i in range(rows-window_size+1):
        for j in range(cols-window_size+1):
          weighted_image[i: i+window_size, j: j+window_size] = diff_array[i,j]

    weighted_images_original.append(weighted_image)

  return weighted_images_original

def find_GA_per_instance_inter(image_info, window_size, sliding_range, weighted_images_original):
  weighted_images_inter = []
  rows, cols = image_info["image_shape"][0],image_info["image_shape"][1]
  for ins in range(1):
    if image_info["instances_categories"][ins] != target_category_id:
      weighted_image = np.zeros((rows, cols))
    else:
      masks_array = np.array(image_info["masks"])
      grad_direction_array = np.array(image_info["grad_fields"]["grad_direction"])
      diff_array = np.zeros((rows-window_size+1, cols-window_size+1))
      for i in range(max(0, sliding_range[1]-window_size+1), min(rows-window_size+1, sliding_range[1]+sliding_range[3])):
        for j in range(max(0, sliding_range[0]-window_size+1), min(cols-window_size+1, sliding_range[0]+sliding_range[2])):
          mask_window = masks_array[:, i:i+window_size, j:j+window_size]
          grad_window = grad_direction_array[:, i:i+window_size, j:j+window_size]
          if np.any(mask_window[ins]):
            diffs = []
            for k in range(image_info["num_instances"]):
              if k != ins and image_info["instances_categories"][k]==target_category_id:
                diff1 = np.std(grad_window[k,:,:] - grad_window[ins,:,:])
                window_1D_k = grad_window[k,:,:].reshape(-1,)
                window_1D_k = np.delete(window_1D_k, np.where(window_1D_k==-4))
                window_1D_i = grad_window[ins,:,:].reshape(-1,)
                window_1D_i = np.delete(window_1D_i, np.where(window_1D_i==-4))
                con_window = np.concatenate((window_1D_k, window_1D_i))
                #diff2 = np.std(con_window)
                #diffs.append(min(diff1, diff2))
                diffs.append(diff1)
            diff_array[i,j] = np.mean(diffs)
      weighted_image = np.zeros((rows,cols))
      for i in range(rows-window_size+1):
        for j in range(cols-window_size+1):
          weighted_image[i: i+window_size, j: j+window_size] = diff_array[i,j]

    weighted_images_inter.append(weighted_image)

  for w in range(len(weighted_images_inter)):
    for k in range(len(weighted_images_original)):
      weighted_images_inter[w] -= weighted_images_original[k]
      weighted_images_inter[w][np.where(weighted_images_inter[w]<0)] = 0
      if k != w and np.logical_and(image_info["masks"][w], image_info["masks"][k]).any() and image_info["instances_categories"][k]==target_category_id and image_info["instances_categories"][w]==target_category_id:
        weighted_images_inter[w][np.logical_and(image_info["masks"][w], image_info["masks"][k])] += overlap_HL
    influence = image_info["grad_fields"]["grad_magnitude"][w]
    influence = influence / np.max(influence)
    weighted_images_inter[w] *= influence
    weighted_images_inter[w] = (weighted_images_inter[w] / (np.max(weighted_images_inter[w])+1e-7))*GA_factor + 1

  return weighted_images_inter

cmap = colors.LinearSegmentedColormap.from_list("mycmap", [(0,0,0,0), 'green', 'blue', 'red', 'yellow'])
def process_dataset(json_data):
  for image in json_data["images"]:
    image_info = {}
    image_info["image_name"] = image["file_name"]
    this_image = mpimg.imread(training_image_path + image_info["image_name"])
    instances_in_this_image = [i for i in json_data['annotations'] if i['image_id'] == image['id']]
    image_info["instances"] = instances_in_this_image
    image_info["instances_ids"] = [i['id'] for i in instances_in_this_image]
    image_info["instances_categories"] = [i["category_id"] for i in instances_in_this_image]
    image_info["num_instances"] = len(instances_in_this_image)
    image_info["image_shape"] = this_image.shape[0:2]
    image_info["masks"] = find_masks(image_info["num_instances"], this_image, instances_in_this_image, image_info["image_shape"])

    image_info["grad_fields"] = find_grad_fields(image_info["masks"], 1)
    weighted_images_original_1 = find_GA_per_instance_original(image_info, window_size, target_category_id)
    weighted_images_inter_1 = find_GA_per_instance_inter(image_info, window_size, target_category_id, weighted_images_original_1)

    image_info["grad_fields"] = find_grad_fields(image_info["masks"], -1)
    weighted_images_original_2 = find_GA_per_instance_original(image_info, window_size, target_category_id)
    weighted_images_inter_2 = find_GA_per_instance_inter(image_info, window_size, target_category_id, weighted_images_original_2)

    weighted_images_inter = [(w1+w2)/2 for w1, w2 in zip(weighted_images_inter_1, weighted_images_inter_2)]

    show_weighted_images_per_instance(weighted_images_inter, this_image)

# json_data = json.load(open(training_json_path))
# process_dataset(json_data)

def find_instances_interaction(json_data):
  # 创建一个字典来存储每个边界框的交集
  intersections = {}
  for image in json_data["images"]:
    instances_in_this_image = [i for i in json_data['annotations'] if i['image_id'] == image['id']]
    this_instance = instances_in_this_image[0]
    # 对于每一个边界框
    for i in instances_in_this_image:
      bbox1 = box(i["bbox"][0] - window_size, i["bbox"][1] - window_size, i["bbox"][0] + i["bbox"][2] + window_size, i["bbox"][1] + i["bbox"][3] + window_size)
      for j in instances_in_this_image:
        if i != j:
          bbox2 = box(j["bbox"][0] - window_size, j["bbox"][1] - window_size, j["bbox"][0] + j["bbox"][2] + window_size, j["bbox"][1] + j["bbox"][3] + window_size)
            # 检查一个边界框是否包含另一个
          if not bbox1.contains(bbox2) and not bbox2.contains(bbox1):
            # 如果没有包含关系，再检查两个边界框是否有交集
            if bbox1.intersects(bbox2):
              # 如果有交集，将它们添加到字典中
              if i['id'] not in intersections:
                intersections[i['id']] = [i]
              intersections[i['id']].append(j)

  # 打印出有交集的边界框
  # for key, values in intersections.items():
  #     print(f"Instance {key} interacts with instances {values}")
  with open('intersections.json', 'w') as f:
    json.dump(intersections, f)
  return intersections

def process_data_with_interaction(json_data, interaction):
  for instance in json_data["annotations"]:
    image_info = {}
    for d in json_data["images"]:
      if d["id"] == instance["image_id"]:
        image_info["image_name"] = d["file_name"]
        break
    #print("image name", image_info["image_name"])
    this_image = mpimg.imread(training_image_path + image_info["image_name"])
    #print("instance id", instance["id"])
    if instance["id"] in interaction:
      instances_with_interaction_in_this_image = interaction[instance["id"]]
      #print("it interacts with instance", interaction[instance["id"]])
      #print(len(instances_with_interaction_in_this_image))
      image_info["instances"] = instances_with_interaction_in_this_image
      image_info["instances_ids"] = [i["id"] for i in instances_with_interaction_in_this_image]
      image_info["num_instances"] = len(instances_with_interaction_in_this_image)
      image_info["instances_categories"] = [i["category_id"] for i in instances_with_interaction_in_this_image]
      #print(image_info["num_instances"])
      image_info["image_shape"] = this_image.shape[0:2]
      image_info["masks"] = find_masks(image_info["num_instances"], this_image, instances_with_interaction_in_this_image, image_info["image_shape"])
      min_x = min(i["bbox"][0] for i in instances_with_interaction_in_this_image)
      min_y = min(i["bbox"][1] for i in instances_with_interaction_in_this_image)
      max_x = max(i["bbox"][0] + i["bbox"][2] for i in instances_with_interaction_in_this_image)
      max_y = max(i["bbox"][1] + i["bbox"][3] for i in instances_with_interaction_in_this_image)
      #print(min_x, min_y, max_x-min_x, max_y-min_y)
      sliding_range = [int(min_x), int(min_y), int(max_x - min_x), int(max_y - min_y)]

      image_info["grad_fields"] = find_grad_fields(image_info["masks"], 1)
      weighted_images_original_1 = find_GA_per_instance_original(image_info, window_size, sliding_range)
      weighted_images_inter_1 = find_GA_per_instance_inter(image_info, window_size, sliding_range, weighted_images_original_1)

      image_info["grad_fields"] = find_grad_fields(image_info["masks"], -1)
      weighted_images_original_2 = find_GA_per_instance_original(image_info, window_size, sliding_range)
      weighted_images_inter_2 = find_GA_per_instance_inter(image_info, window_size, sliding_range, weighted_images_original_2)

      weighted_images_inter = [(w1+w2)/2 for w1, w2 in zip(weighted_images_inter_1, weighted_images_inter_2)]
    else:
      weighted_images_inter = [np.zeros(this_image.shape[0:2])]
    #print("instance_id", instance["id"])

    #show_weighted_images_per_instance(weighted_images_inter, this_image)
    save_weighted_npy_per_instance(weighted_images_inter, weighted_image_path, [instance["id"]])


json_data = json.load(open(training_json_path))
interaction = find_instances_interaction(json_data)
#with open('intersections.json', 'r') as f:
#    intersections = json.load(f)
process_data_with_interaction(json_data, interaction)