python-3.x - 力矩保持阈值 : Trilevel Case

Question

我目前正在实现 Tsai 提出的分割方法对于 Python 3.6 中 8 位图像中的三级情况。这是代码，考虑到论文附件中解释的数学关系:

 ### EXAMPLE GIVEN BY TSAI ### 
 
 # This is the "dummy" example described in the paper
 img = np.array([[10, 8, 10, 9, 20, 21,32,30,40,41,41,40],
                   [12, 10, 11, 10, 19, 20, 30, 28, 38, 40, 40, 39],
                   [10, 9, 10, 8, 20, 21, 30, 29, 42, 40, 40, 39],
                   [11, 10, 9, 11, 19, 21, 31, 30, 40, 42,38, 40]])
    
    bins_hist = list(range(0,257)) 
    histogram = np.zeros((256,4))  # We generate the counts for each grey value, 3rd column probabilities and 4th column cumulative sum
    histogram[:,0] = np.arange(256)  # First column contains the intensity levels
    
    hist_i = np.histogram(img, bins=bins_hist)
    counts = hist_i[0]
    histogram[:,1] = np.add(histogram[:,1], counts)
   
    row, col = img.shape
    total_voxels = row * col
    
    histogram[:,2] = histogram[:,1]/total_voxels  # Relative frequency at each GV
    histogram[:,3] = np.cumsum(histogram[:,2])  # Cumulative of the relative frequency


 ###  THRESHOLDS CALCULATION BEGIN ###
    
    # 0 Moment = 1
    m0 = 1.0
    
    # 1st Moment 
    m1 = np.cumsum((histogram[:,0])*histogram[:,2])[-1]    
    
    # 2nd Moment
    m2 = np.cumsum((histogram[:,0]**2)*histogram[:,2])[-1]  # Take the last value
    
    # 3rd Moment
    m3 = np.cumsum((histogram[:,0]**3)*histogram[:,2])[-1]
    
    # 4th Moment
    m4 = np.cumsum((histogram[:,0]**4)*histogram[:,2])[-1] 
    
    # 5th Moment
    m5 = np.cumsum((histogram[:,0]**5)*histogram[:,2])[-1]
    
    # Now we must find the value in the binary image that preserves these moments
    
    
    # We solve the equalities --> For solutions refer to Paper Annex A.2
    cd = (m0*m2*m4) + (m1*m3*m2) + (m1*m3*m2) - (m2*m2*m2) - (m1*m1*m4) - (m3*m3*m0)
    c0 = ((-m3*m2*m4) + (-m4*m3*m2) + (m1*m3*-m5) - (-m5*m2*m2) - (-m4*m1*m4) - (m3*m3*-m3)) / cd
    c1 = ((m0*-m4*m4) + (m1*-m5*m2) + (-m3*m3*m2) - (m2*-m4*m2) - (m1*-m3*m4) - (-m5*m3*m0)) / cd
    c2 = ((m0*m2*-m5) + (m1*m3*-m3) + (m1*-m4*m2) - (m2*m2*-m3) - (m1*m1*-m5) - (m3*-m4*m0)) /cd
    
    a1 = c0/2 - c1*c2/6 + (c2**3)/27
    a2 = (c0/2 - c1*c2/6 + (c2**3)/27)**2
    a3 = (c1/3 - (c2**2)/9)**3
    
    a = (a1 - cmath.sqrt(a2 + a3))**1/3 
    
    b = -(c1/3 - (c2**2)/9)/a
    w1 = -0.5 + 1j * (math.sqrt(3)/2)
    w2 = -0.5 - 1j * (math.sqrt(3)/2)
    
    z0 = -c2/3 - a - b
    z1 = -c2/3 - w1*a - w2*b
    z2 = -c2/3 - w2*a - w1*b
   
    pd = (z1*z2**2) + (z2*z0**2) + (z0*z1**2) - (z0**2*z1) - (z0*z2**2) - (z1**2*z2)
    p0 = ((m0*z1*z2**2) + (m1*z1**2) + (z2*m2) - (m2*z1) - (m1*z2**2) - (z1**2*z2*m0)) /pd
    p1 = ((m1*z2**2) + (z0*m2) + (m0*z2*z0**2) - (z0**2*m1) - (z0*m0*z2**2) - (m2*z2)) / pd # Fraction of the below-threshold pixels in the binary histogram

    th1 = 0.0  # First threshold in the trimodal histogram
    th2 = 0.0  # Second threshold 
    dist1 = 10000000
    dist2 = 10000000
    
    for i in range(254):
        for j in range(i+1, 255):
            # Select threshold --> closest to p0 from the normlaized histogram
            p0_orig = histogram[i,3]  # Take the cumulative relative frequency at the value p0
            p1_orig = histogram[j, 3] - histogram[i,3]
            dist_i = abs(p0 - p0_orig)
            dist_j = abs(p1 - p1_orig)
            
            if dist_i < dist1 and dist_j < dist2:  # This one was the one mentioned by Tsai ("Minimize the distance")
                
                print(i,j,dist_i, dist_j)
                dist1 = dist_i
                dist2 = dist_j
                th1 = i 
                th2 = j

然而，当我考虑到图 1 中描述的“虚拟”图像而应用它来重现他的结果时，我没有获得相同的阈值(在我的情况下，我得到 12 和 29，而不是 18 和 30)。
有没有人有这种方法的经验，可以帮助我找出我的代码有什么问题？

Answer 1

可能有一些错别字，例如:a = (a1 - cmath.sqrt(a2 + a3)) ** 1 / 3这实际上应该是:a = (a1 - cmath.sqrt(a2 + a3)) ** (1 / 3)并建议使用numpy线性代数的函数来避免一些错误，比如npmpy.dot(x,y), numpy.vdot(x,y)等等。这也使代码更具可读性。