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?
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.