语义分割后处理，单细胞分析和目标追踪

用于细胞分割之后的随时间变化的单细胞追踪，如果 optic flow就不够精确，所以我做了这个，基本原理非常弱智，可以自定义的threshhold很多，而且可以很快划掉比较快。后面我加了匈牙利算法，请参考下一篇关于匈牙利算法的文章

import cv2 as cv
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure, color
import math
from sklearn import linear_model  # 表示，可以调用sklearn中的linear_model模块进行线性回归。
import tifffile
import os
from PIL import Image
from glob import glob
from tqdm import tqdm
import pandas as pd
import seaborn as sns

import time
import skimage

"""
工具 二值化
thresh=0 除背景以外为1
thresh=1 border + cyto
thresh=2 cyto
thresh=3 nuclei
显示单个细胞
"""


def get_img_binary(pred, thresh):
    ret, thresh_pred_1 = cv.threshold(pred, thresh, 1, cv.THRESH_BINARY)
    return thresh_pred_1


def represent_single_cell_from(pred_info, num):
    if num is not None:
        coords = pred_info[num][5]
        ohhh_img = np.zeros((1200, 1200))
        for i in coords:
            x = int(i[0])
            y = int(i[1])
            ohhh_img[x, y] = 1
        plt.imshow(ohhh_img)
    else:
        print("Alert! : no matching objects,num is None!")


"""
筛选细胞面积，
闭开运算分割细胞去除噪音

获得筛选后标签图片
在这种情况下，我们得到更多的单细胞，但是损失了细节
    开运算：先腐蚀后膨胀, 去除噪声，去除白色小点、空洞
    闭运算：先膨胀后腐蚀, 用来填充前景物体的小黑点
    形态学梯度：膨胀减去腐蚀, 可以得到前景物体的轮廓
    礼帽：原图减去开运算
    黑帽：闭运算减去原图
"""


def get_labeled_pred(pred_img, thresh=1):
    ret, thresh_pred = cv.threshold(np.array(pred_img), thresh, 1, cv.THRESH_BINARY)
    # 形态学调整
    # open 去噪
    # close 闭口
    kernel = np.ones((3, 3), np.uint8)
    closing = cv.morphologyEx(thresh_pred, cv.MORPH_CLOSE, kernel)
    kernel2 = np.ones((8, 8), np.uint8)
    opening = cv.morphologyEx(closing, cv.MORPH_OPEN, kernel2)

    labels_objects = measure.label(opening, background=0, connectivity=2)

    properties = measure.regionprops(labels_objects)
    object_areas = []
    object_labels = []
    object_coords = []
    for prop in properties:
        object_areas.append(prop.area)
        object_labels.append(prop.label)
        object_coords.append(prop.coords)
    valid_labels = []
    for label1, area in enumerate(object_areas):
        # ============================
        # 确定筛选范围，需要随时调整的
        # ============================
        if area < 2800 and area > 600:
            valid_labels.append(label1)
    valid_coords = []
    for label2, coord in enumerate(object_coords):
        for label1 in valid_labels:
            if label2 == label1:
                valid_coords.append(object_coords[label2])
    imgnew = np.zeros((1200, 1200))
    for k in valid_coords:
        for i in k:
            imgnew[i[0], i[1]] = 1
    valid_labels_objects = measure.label(imgnew, background=0, connectivity=2)
    return valid_labels_objects


"""
在label删选过的基础上获取信息
获取的信息注意cell——0对应标签图片的label——1
0 label
1 面积 area
2 周长 perimeter
3 重心 centroid
4 bbox 
5 坐标 coords
6 离心 en
"""


def get_labeled_info(labeled_pred):
    properties = measure.regionprops(labeled_pred)
    object_areas = []
    object_perimeters = []  # 周长
    object_centroids = []  # 重心点
    object_bboxs = []
    object_coords = []  # 区域内像素点坐标
    object_labels = []
    object_eccentricities = []
    object_images = []
    valid_para_of_1image = []
    for prop in properties:
        object_areas.append(prop.area)
        object_perimeters.append(prop.perimeter)
        object_centroids.append(prop.centroid)
        object_bboxs.append(prop.bbox)
        object_coords.append(prop.coords)
        object_labels.append(prop.label)
        object_eccentricities.append(prop.eccentricity)
        object_images.append(prop.image)
    ziped_para_of_image = list(zip(object_labels,  # 0
                                   object_areas,  # 1 面积
                                   object_perimeters,  # 2周长
                                   object_centroids,  # 3重心点
                                   object_bboxs,
                                   object_coords,  # 5区域内像素点坐标
                                   object_eccentricities,
                                   object_images))
    return ziped_para_of_image


"""
获得筛选过的labeled pred list
同步获得他们的图片信息
"""


def get_labeled_preds_and_info(pred_list):
    labeled_preds = []
    preds_info = []
    for pred in tqdm(pred_list):
        labeled_pred_cell = get_labeled_pred(pred, thresh=1)
        labeled_preds.append(labeled_pred_cell)
        preds_info.append(get_labeled_info(labeled_pred_cell))
    print(len(labeled_preds), len(preds_info))
    return labeled_preds, preds_info


def make_small_img(target_cell_num, labeled_pred_img, pred_info, size):
    x = int(pred_info[target_cell_num][3][0])
    y = int(pred_info[target_cell_num][3][1])
    small_img = labeled_pred_img[np.max([0, x - size]):np.min([x + size, 1200]),
                np.max([0, y - size]):np.min([y + size, 1200])]

    return small_img


def get_first_pred(labeled_preds):
    pred_1 = labeled_preds[0]
    plt.figure(figsize=(24, 24))
    plt.imshow(pred_1)
    labeled_pred1 = get_labeled_pred(labeled_preds[0], 1)
    print(np.max(labeled_pred1))
    pred1_info = get_labeled_info(labeled_pred1)
    print(len(pred1_info))
    plt.figure(figsize=(24, 24))
    plt.imshow(labeled_pred1)
    return labeled_pred1, pred1_info


def find_color_label(num, labeled_pred1, labeled_pred2, pred1_info, pred2_info, size):
    # 以第一张图的重心为重心，划出目标图局部区域
    simg2 = make_small_img(num, labeled_pred2, pred1_info, size)
    # plt.imshow(simg2)
    cell_num = np.unique(simg2)[1:] - 1

    distances = []
    for label in cell_num:
        x = int(pred2_info[label][3][0])
        y = int(pred2_info[label][3][1])
        distance = math.sqrt(((int(pred1_info[num][3][0]) - x) ** 2 + (int(pred1_info[num][3][1]) - y) ** 2))

        distances.append(distance)
    # print(distances)
    target = None
    if np.size(distances) == 0:
        target = None
    else:
        if np.min(distances) < 20:
            target = cell_num[np.argmin(distances)]

    return target


def get_same_cell_list(num, labeled_pred1, labeled_preds, pred1_info, preds_info, size):
    same_cell_list = []
    for i, img in enumerate(labeled_preds):
        same_cell = find_color_label(num, labeled_pred1, img, pred1_info, preds_info[i], size)
        same_cell_list.append(same_cell)
    return same_cell_list


"""
用同一套prediction计算两个相应的FRET
相减 获得细胞亮度
"""


def calculate_brightness(num, target_imgfret1, target_imgfret2, one_pred_info, one_labeled_pred):
    box = one_pred_info[num][4]

    onecellmask = one_pred_info[num][7]
    onecellarea = one_pred_info[num][1]

    small_target_pred = one_labeled_pred[box[0]:box[2], box[1]:box[3]]
    maskofallcell = get_img_binary(np.float64(small_target_pred), 0)
    bg_reversmask = [[1 - i for i in j] for j in maskofallcell]
    bg_area = np.sum(bg_reversmask)

    boxinfret1 = target_imgfret1[box[0]:box[2], box[1]:box[3]]
    validfret1 = boxinfret1 * onecellmask
    outsiderfret1 = boxinfret1 * bg_reversmask

    boxinfret2 = target_imgfret2[box[0]:box[2], box[1]:box[3]]
    validfret2 = boxinfret2 * onecellmask
    outsiderfret2 = boxinfret2 * bg_reversmask
    f1 = np.sum(validfret1) / onecellarea - np.sum(outsiderfret1) / bg_area
    f2 = np.sum(validfret2) / onecellarea - np.sum(outsiderfret2) / bg_area
    return (f2 / f1)


def calculate_brightness_oldversion(num, target_imgfret1, target_imgfret2, One_pred_info):
    bright1 = []
    bright2 = []
    pppsjfk = target_imgfret1.copy()
    pppsjfk2 = target_imgfret2.copy()
    for i in One_pred_info[num][5]:
        x = i[0]
        y = i[1]
        bright1.append(pppsjfk[x, y])
        bright2.append(pppsjfk2[x, y])
    f1 = sum(bright1) / One_pred_info[num][1]
    f2 = sum(bright2) / One_pred_info[num][1]
    return (f1 / f2)


def get_brightness_area_perimeter_eccentricities(num, labeled_pred1, labeled_preds, pred1_info, preds_info, size,
                                                 image_FRET1, image_FRET2, add_agent_frame, figure):
    same_cell_list = get_same_cell_list(num, labeled_pred1, labeled_preds, pred1_info, preds_info, size)
    sequence_of_imgs = range(len(pred_FRET1))
    how_many_valida_cells = sum(x is not None for x in same_cell_list)
    zipped = list(zip(sequence_of_imgs, same_cell_list))
    fretlist = []
    fretlist_after_agent = []
    arealist = []
    arealist_after_agent = []
    perimeterlist = []
    perimeterlist_after_agent = []
    eccentricitieslist = []
    eccentricitieslist_after_agent = []
    if how_many_valida_cells > 0.8 * len(pred_FRET1):
        for ll, (v, k) in enumerate(zipped):
            if k is not None:
                brightnessRatio = calculate_brightness(k, image_FRET1[v], image_FRET2[v], preds_info[v],
                                                       labeled_preds[v])
                fretlist.append(brightnessRatio)
                area = preds_info[v][k][1]
                arealist.append(area)
                perimeter = preds_info[v][k][2]
                perimeterlist.append(perimeter)
                eccentricity = preds_info[v][k][6]
                eccentricitieslist.append(eccentricity)
                if ll > add_agent_frame and ll < add_agent_frame + 40:
                    fretlist_after_agent.append(brightnessRatio)
                    area = preds_info[v][k][1]
                    arealist_after_agent.append(area)
                    perimeter = preds_info[v][k][2]
                    perimeterlist_after_agent.append(perimeter)
                    eccentricity = preds_info[v][k][6]
                    eccentricitieslist_after_agent.append(eccentricity)

        if figure:
            x = range(how_many_valida_cells)
            y1 = np.array(fretlist)
            y2 = np.array(arealist)

            y3 = np.array(perimeterlist)
            y4 = np.array(eccentricitieslist)
            plt.subplot(4, 1, 1)
            plt.title("FRET")
            plt.plot(x, y1)
            plt.subplot(4, 1, 2)
            plt.title("AREA")
            plt.plot(x, y2)
            plt.subplot(4, 1, 3)
            plt.title("Perimeter")
            plt.plot(x, y3)
            plt.subplot(4, 1, 4)
            plt.title("eccentricities")
            plt.plot(x, y4)

        return fretlist_after_agent, arealist_after_agent, perimeterlist_after_agent, eccentricitieslist_after_agent, fretlist, arealist, perimeterlist, eccentricitieslist
    else:
        return None


def get_num_and_slops_of_cells_in_first_image(labeled_pred1, labeled_preds, pred1_info, preds_info, image_FRET1,
                                              image_FRET2, add_agent_frame):
    slops = []
    for i in tqdm(range(len(preds_info[0]))):
        # for i in tqdm(range(least_cell_numbers(labeled_preds))):
        M = get_brightness_area_perimeter_eccentricities(i, labeled_pred1, labeled_preds, pred1_info, preds_info, 80,
                                                         image_FRET1, image_FRET2, add_agent_frame, figure=False)

        if M is not None:
            a, b, c, d, a1, b1, c1, d1 = M
            model = linear_model.LinearRegression()
            x = range(len(a))
            x = np.reshape(x, (-1, 1))
            model.fit(x, np.array(a))
            slop = model.coef_
            slops.append(slop)
    return len(preds_info[0]), slops


def check_same_cell_centro_and_position_in_two_preds(num_of_cell, num_of_img, num_of_img2, preds_info):
    print(preds_info[num_of_img][num_of_cell][3])
    neww = np.zeros((1200, 1200))
    for i in preds_info[num_of_img][num_of_cell][5]:
        neww[i[0], i[1]] = 1
    plt.figure(figsize=(24, 24))
    plt.subplot(3, 1, 1)
    plt.imshow(neww)
    print(preds_info[num_of_img2][num_of_cell][3])
    neww3 = np.zeros((1200, 1200))
    for i in preds_info[num_of_img2][num_of_cell][5]:
        neww3[i[0], i[1]] = 1
    plt.subplot(3, 1, 2)
    plt.imshow(neww3)


def make_box_fig(slops1, slops2):
    plt.style.use('_mpl-gallery')

    # make data:
    np.random.seed(10)
    D = slops1, slops2

    # plot
    fig, ax = plt.subplots()
    VP = ax.boxplot(D, positions=[2, 4], widths=1.5, patch_artist=True,
                    showmeans=False, showfliers=False,
                    medianprops={"color": "white", "linewidth": 0.5},
                    boxprops={"facecolor": "C0", "edgecolor": "white",
                              "linewidth": 0.5},
                    whiskerprops={"color": "C0", "linewidth": 1.5},
                    capprops={"color": "C0", "linewidth": 1.5})

    # ax.set(xlim=(0, 6), xticks=np.arange(1, 6),
    #        ylim=(0, 6), yticks=np.arange(1, 6))

    plt.show()


def dataloader_pred_fret1_fret2(root_dir, seg_folder_name, extend1, extend2):

    pred_FRET1_path = glob(os.path.join(root_dir, seg_folder_name, '*' + extend1 + ".png"))
    FRET1_path = glob(os.path.join(root_dir, '*' + extend1 + ".tif"))
    FRET2_path = glob(os.path.join(root_dir, '*' + extend2 + ".tif"))
    pred_FRET1 = [Image.open(img) for img in pred_FRET1_path]
    image_FRET1 = [tifffile.imread(img) for img in FRET1_path]
    image_FRET2 = [tifffile.imread(img) for img in FRET2_path]
    print(
        f'预测fret1图片数：{len(pred_FRET1)},,原始fret1图片数：{len(FRET1_path)},原始fret2图片数：{len(FRET2_path)},', )
    return pred_FRET1, image_FRET1, image_FRET2


"""
追踪同一细胞的可视化
visuali the tracked cell
"""


def visualization_of_tracking_sys(test_cell_num, show_it_on_which_image, labeled_preds, preds_info):
    # print(len(pred1_info),len(pred2_info),len(pred3_info))
    # 要查看的细胞
    iii = show_it_on_which_image  # 和第一张图片比较的图片的序号
    labeled_pred1 = labeled_preds[0]
    pred1_info = preds_info[0]
    target_in_pred2 = find_color_label(test_cell_num, labeled_pred1, labeled_preds[iii], pred1_info, preds_info[iii],
                                       size=80)
    plt.figure(figsize=(18, 18))
    plt.subplot(1, 2, 1)
    represent_single_cell_from(pred1_info, test_cell_num)
    plt.subplot(1, 2, 2)
    represent_single_cell_from(preds_info[iii], target_in_pred2)


"""
画出箱线图
draw box figure
"""


def get_box_plot(slops1, slops2):
    len1 = len(np.array(slops1))
    len2 = len(np.array(slops2))
    list1 = np.ones((len1,))
    list2 = np.ones((len2,)) * 2
    num_list = pd.concat((pd.DataFrame(list1), pd.DataFrame(list2)), axis=0)
    data = pd.concat([pd.DataFrame(slops1), pd.DataFrame(slops2)], axis=0)
    num_list = num_list.reset_index()
    data = data.reset_index()
    data_all = pd.concat([pd.DataFrame(num_list), pd.DataFrame(data)], axis=1)
    data_all.columns = (["index1", "exp", "index2", "slopes"])
    # print(data_all.head())
    sns.boxenplot(x="exp", y="slopes", data=data_all)