如何计算 1D 中的 Otsu 阈值答案

【问题标题】：how to calculate Otsu threshold in 1D如何计算 1D 中的 Otsu 阈值
【发布时间】：2019-03-14 19:54:58
【问题描述】：

我正在尝试识别分析化学数据中的双峰分布。每个数据集是来自 GC-MS 的特定化合物的 3~70 个保留时间列表。某些化合物的 RT 呈双峰分布，其中库搜索已将相同的身份分配给具有不同 RT 的数据中的两个或多个不同特征。这对于具有非常相似质谱的异构体和其他化合物对来说是很常见的。例如。这是显示双峰分布的一种化合物的 RT 直方图。

我想计算 Otsu 阈值以尝试定义双峰数据（也有多峰分布，但一次只有一步）。我很难理解the Wikipedia article 的计算，但文本表明可以通过找到最小的类内方差来找到阈值。因此，我尝试从 RT 列表中进行计算，如下所示：

a = list(d['Component RT'])
n = len(a)
b = [a.pop(0)]

varA = []
varB = []

for i in range(1,n-2):
    b.append(a.pop(0))
    varA.append(statistics.stdev(a)**2)
    varB.append(statistics.stdev(b)**2)

我是否认为如果我绘制上述数据的方差总和，我应该能够将 Otsu 阈值识别为最小值？

在此示例中，阈值很明显，大约有 35 个值可供使用。对于大多数化合物，值较少（通常

【问题讨论】：

标签： threshold

【解决方案1】：

从这个answer我们可以推断出一维情况如下：

import numpy as np

# source: https://stackoverflow.com/a/50796152

def otsu1d(array1d, maxoutval = 255):
    pixel_number = array1d.shape[0]
    mean_weight = 1.0/pixel_number
    his, bins = np.histogram(array1d, np.arange(0,np.max(array1d)+1))
    final_thresh = -1
    final_value = -1
    intensity_arr = np.arange(np.max(array1d))
    for t in bins[1:-1]: # This goes from 1 to 254 uint8 range (Pretty sure wont be those values)
        pcb = np.sum(his[:t])
        pcf = np.sum(his[t:])
        Wb = pcb * mean_weight
        Wf = pcf * mean_weight

        mub = np.sum(intensity_arr[:t]*his[:t]) / float(pcb)
        muf = np.sum(intensity_arr[t:]*his[t:]) / float(pcf)
        #print mub, muf
        value = Wb * Wf * (mub - muf) ** 2

        if value > final_value:
            final_thresh = t
            final_value = value
    final_arr = array1d.copy()
    print(final_thresh)
    final_arr[array1d > final_thresh] = maxoutval
    final_arr[array1d < final_thresh] = 0
    return final_arr

【讨论】：