Moment Preserving Thresholding: Trilevel Case

Question

I'm currently implementing the segmentation method proposed by Tsai for the trilevel case in an 8-bit image in Python 3.6. This is the code, considering the mathematical relations explained in the paper Annex:


 ### 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

However, when I apply it to reproduce his results considering a "dummy" image described in Figure 1, I do not obtain the same threshold (in my case I get 12 and 29, instead of 18 and 30).

Does anybody have experience with this method and can help me figure out what's the problem with my code?

Hi! I have just updated the question by including the array that Tsai provides in his paper. — pedro_galher, Sep 21 '21 at 15:35
The solution is likely not unique. Maybe your solution and the author’s are equivalent, or produce the same error. And of course it is also possible that there’s a typo in the paper, or that the authors didn’t implement their own math right. Not saying this is likely the case, but I’ve seen it happen. — Cris Luengo, Sep 24 '21 at 22:00
@CrisLuengo yes, you are completely right. Actually, you gave me a good idea and I will double-check that maybe my result and Tsai's one produce the same moments (up to 5th degree). Thanks! — pedro_galher, Sep 24 '21 at 22:08

score 1 · Accepted Answer · answered Sep 25 '21 at 06:50

1

There may be some typos like: a = (a1 - cmath.sqrt(a2 + a3)) ** 1 / 3

which actually should be: a = (a1 - cmath.sqrt(a2 + a3)) ** (1 / 3)

And suggest to use the functions of numpy linear algebraic to avoid some mistakes, such as npmpy.dot(x,y), numpy.vdot(x,y) and so on. which also make the code much more readable.

answered Sep 25 '21 at 06:50

NicoNing

3,076
12
23

Thanks for your answer. Actually in this code no algebraic functions are required as there is no matrix multiplication rather element-wise multiplication. I'll revise the code again to try and find more typos. Thanks! – pedro_galher Sep 25 '21 at 09:40
1

I forgot to add that you spotted a very important typo! Thanks to that the thresholds I get (19 and 30) are very close to the Tsai ones (18 and 30) – pedro_galher Sep 25 '21 at 09:52
After correcting the code following your comment and comparing the rest of the parameters (p0, p1, p1, z0, z1 and z2) with what Tsai obtains on his paper, I think that this difference of just one pixel is due to rounding errors as Tsai gives everything without decimals but I did not do any rounding in the code. I believe that correcting the typo made the code work. Thank you again! – pedro_galher Sep 25 '21 at 21:06

Moment Preserving Thresholding: Trilevel Case

1 Answers1