用于腺体分割的标准参数：对象级Dice、对象级豪斯多夫距离

在网上都没找到ObjectDice、ObjectHausdorffDistance、F1score的代码，在这里给大家贴上

import numpy as np
from scipy import stats
from sklearn.neighbors import NearestNeighbors


def ObjectHausdorff(S=None, G=None):
    S = np.array(S).astype(np.uint8)
    G = np.array(G).astype(np.uint8)

    totalAreaS = (S > 0).sum()
    totalAreaG = (G > 0).sum()
    listLabelS = np.unique(S)
    listLabelS = np.delete(listLabelS, np.where(listLabelS == 0))
    listLabelG = np.unique(G)
    listLabelG = np.delete(listLabelG, np.where(listLabelG == 0))

    temp1 = 0
    for iLabelS in range(len(listLabelS)):
        Si = (S == listLabelS[iLabelS])
        intersectlist = G[Si]
        if intersectlist.any():
            indexGi = stats.mode(intersectlist).mode
            Gi = (G == indexGi)
        else:
            tempDist = np.zeros((len(listLabelG), 1))
            for iLabelG in range(len(listLabelG)):
                Gi = (G == listLabelG[iLabelG])
                tempDist[iLabelG] = Hausdorff(Gi, Si)
            minIdx = np.argmin(tempDist)
            Gi = (G == listLabelG[minIdx])
        omegai = Si.sum() / totalAreaS
        temp1 = temp1 + omegai * Hausdorff(Gi, Si)

    temp2 = 0
    for iLabelG in range(len(listLabelG)):
        tildeGi = (G == listLabelG[iLabelG])
        intersectlist = S[tildeGi]
        if intersectlist.any():
            indextildeSi = stats.mode(intersectlist).mode
            tildeSi = (S == indextildeSi)
        else:
            tempDist = np.zeros((len(listLabelS), 1))
            for iLabelS in range(len(listLabelS)):
                tildeSi = (S == listLabelS[iLabelS])
                tempDist[iLabelS] = Hausdorff(tildeGi, tildeSi)
            minIdx = np.argmin(tempDist)
            tildeSi = (S == listLabelS[minIdx])
        tildeOmegai = tildeGi.sum() / totalAreaG
        temp2 = temp2 + tildeOmegai * Hausdorff(tildeGi, tildeSi)

    objHausdorff = (temp1 + temp2) / 2
    return objHausdorff


def Hausdorff(S=None, G=None, *args, **kwargs):
    S = np.array(S).astype(np.uint8)
    G = np.array(G).astype(np.uint8)

    listS = np.unique(S)
    listS = np.delete(listS, np.where(listS == 0))
    listG = np.unique(G)
    listG = np.delete(listG, np.where(listG == 0))

    numS = len(listS)
    numG = len(listG)
    if numS == 0 and numG == 0:
        hausdorffDistance = 0
        return hausdorffDistance
    else:
        if numS == 0 or numG == 0:
            hausdorffDistance = np.Inf
            return hausdorffDistance

    y = np.where(S > 0)
    x = np.where(G > 0)

    x = np.vstack((x[0], x[1])).transpose()
    y = np.vstack((y[0], y[1])).transpose()

    nbrs = NearestNeighbors(n_neighbors=1).fit(x)
    distances, indices = nbrs.kneighbors(y)
    dist1 = np.max(distances)

    nbrs = NearestNeighbors(n_neighbors=1).fit(y)
    distances, indices = nbrs.kneighbors(x)
    dist2 = np.max(distances)

    hausdorffDistance = np.max((dist1, dist2))
    return hausdorffDistance


def F1score(S=None, G=None):
    S = np.array(S).astype(np.uint8)
    G = np.array(G).astype(np.uint8)

    unique_values_S = np.unique(S)
    unique_values_S = np.delete(unique_values_S, np.where(unique_values_S == 0))

    unique_values_G = np.unique(G)
    unique_values_G = np.delete(unique_values_G, np.where(unique_values_G == 0))

    precision_list = []
    recall_list = []

    for value_S in unique_values_S:
        for value_G in unique_values_G:
            if value_S != 0:
                SegObj = (S == value_S)
                GTObj = (G == value_G)

                overlap = np.logical_and(SegObj, GTObj)
                areaOverlap = overlap.sum()
                areaGTObj = GTObj.sum()

                if areaGTObj > 0:
                    precision = areaOverlap / areaGTObj
                    recall = areaOverlap / GTObj.sum()

                    precision_list.append(precision)
                    recall_list.append(recall)

    precision = np.sum(precision_list) / len(precision_list) if precision_list else 0
    recall = np.sum(recall_list) / len(recall_list) if recall_list else 0

    # Handle the case when precision + recall is zero
    if precision + recall == 0:
        return 0

    f1_score = (2 * precision * recall) / (precision + recall)

    return f1_score


def ObjectDice(S, G):
    S = np.array(S).astype(np.uint8)
    G = np.array(G).astype(np.uint8)

    totalAreaG = (G > 0).sum()
    listLabelS = np.unique(S)
    listLabelS = np.delete(listLabelS, np.where(listLabelS == 0))
    numS = len(listLabelS)
    listLabelG = np.unique(G)
    listLabelG = np.delete(listLabelG, np.where(listLabelG == 0))
    numG = len(listLabelG)

    if numS == 0 and numG == 0:
        return 1
    elif numS == 0 or numG == 0:
        return 0

    temp1 = 0
    totalAreaS = (S > 0).sum()
    for iLabelS in range(len(listLabelS)):
        Si = (S == listLabelS[iLabelS])
        intersectlist = G[Si]
        if intersectlist.any():
            indexG1 = stats.mode(intersectlist).mode
            Gi = (G == indexG1)
        else:
            Gi = np.zeros(G.shape)

        omegai = Si.sum() / totalAreaS
        temp1 += omegai * Dice(Gi, Si)

    temp2 = 0
    totalAreaG = (G > 0).sum()
    for iLabelG in range(len(listLabelG)):
        tildeGi = (G == listLabelG[iLabelG])
        intersectlist = S[tildeGi]
        if intersectlist.any():
            indextildeSi = stats.mode(intersectlist).mode
            tildeSi = (S == indextildeSi)  # np logical and?
        else:
            tildeSi = np.zeros(S.shape)

        tildeOmegai = tildeGi.sum() / totalAreaG
        temp2 += tildeOmegai * Dice(tildeGi, tildeSi)

    return (temp1 + temp2) / 2


def Dice(A, B):
    intersection = np.logical_and(A, B)
    return 2. * intersection.sum() / (A.sum() + B.sum())