【发布时间】：2021-07-29 05:31:01
【问题描述】：
我目前正在实现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）。
有没有人有这种方法的经验并且可以帮助我找出我的代码有什么问题？
【问题讨论】：
-
你能把虚拟图像贴在这里吗？
-
嗨！我刚刚通过包含 Tsai 在他的论文中提供的数组来更新问题。
-
解决方案可能不是唯一的。也许您的解决方案和作者的解决方案是等效的，或者产生相同的错误。当然，也有可能是论文中有错别字，或者作者没有正确实现自己的数学。并不是说这很可能是这种情况，但我已经看到它发生了。
-
@CrisLuengo 是的，你完全正确。实际上，你给了我一个好主意，我会仔细检查我的结果和 Tsai 的结果是否产生相同的时刻（最高 5 级）。谢谢！
标签： python-3.x image-processing image-segmentation image-thresholding